REFMAP Listening test 1 Parts A & B analysis: Random forest variable importance identification - Revision 1

This notebook analyses both parts of the data in terms of variable importance, using a random forest model based on conditional inference trees and a conditional permutation variable importance algorithm.

Setup

Load packages

# load packages
library(tidyr)
library(ggplot2)
library(party)
library(conflicted)
library(tidyverse)
library(openxlsx)
library(caret)
library(viridis)
library(cowplot)
library(permimp)
# set package parameters
theme_set(theme_bw())

# plot colour scheme

mycolourlist = list(c(0, 102, 255), c(0, 204, 153), c(255, 0, 102), c(74, 111, 152), c(251, 164, 49), c(204, 153, 255), c(90, 192, 255), c(80, 245, 233), c(255, 90, 192), c(164, 201, 242), c(255, 254, 139), c(255, 243, 255))
mycolours = matrix()

for (ii in 1:length(mycolourlist)){
  mycolours[ii] = rgb(mycolourlist[[ii]][1]/255,
                      mycolourlist[[ii]][2]/255,
                      mycolourlist[[ii]][3]/255)
}

# toggle to save plots
saveplots = FALSE

if (saveplots){
  # set output plot directory
  choose.files(caption="Just cancel this", filters=matrix(data=c(" ", " "), ncol=2))  # workaround for bug in RTerm choose.dir
  outFigPath <- utils::choose.dir(caption="Select output folder to save plots '03 Experiment\\Experiment 1\\Analysis\\Plots'")
  
  if (!dir.exists(file.path(outFigPath, "svg"))){dir.create(file.path(outFigPath, "svg"))}
  if (!dir.exists(file.path(outFigPath, "pdf"))){dir.create(file.path(outFigPath, "pdf"))}
  
}

# toggle to save data
savedata = TRUE

if (savedata){
  # set output plot directory
  if (saveplots==FALSE){
    choose.files(caption="Just cancel this", filters=matrix(data=c(" ", " "), ncol=2))  # workaround for bug in RTerm choose.dir
  }
  outDataPath <- utils::choose.dir(caption="Select output folder to save data '03 Experiment\\Experiment 1\\Analysis\\R'")
}
 

Import data and wrangle

stimDatapath <- utils::choose.files(caption=r"(Select refmap_listest1_testdata_ByStim.csv from 03 Experiment\Experiment 1\Analysis\PostProcess)",
                                     filters=matrix(data=c("refmap_listest1_testdata_ByStim.csv", "refmap_listest1_testdata_ByStim.csv"), ncol=2))

stimData <- utils::read.csv(stimDatapath, header=TRUE)

colnames(stimData)[1] <- "Stimulus"

# make response proportions into percentages
stimData[['HighAnnoyPc']] <- stimData[['HighAnnoyProp']]*100
stimData[['dHighAnnoyPc']] <- stimData[['dHighAnnoyProp']]*100

# function to encode categorical to ordinal numeric variables
encode_ordinal <- function(x, order=unique(x)) {
  x <- as.numeric(factor(x, levels=order, exclude=NULL, order=TRUE))
  x
}

# definition of ordinal variable levels
SNRCats <- c("No UAS", "-16", "-10", "-4", "2", "8")
UASLAeqCats <- c("No UAS", "42", "48", "54", "60")

The aggregated data by stimulus are assigned to a dataframe, relevant categorical variables are converted to ordinal, and then the variable subset of interest is selected, NA rows dropped (ie, the ‘no UAS’ stimuli, as the conditional variable importance algorithm cannot currently handle NA values, which are present in all the UAS dB metrics), and a formula assigned.

stimDataNum <- data.frame()

stimDataNum <- cbind(stimData[, 'Stimulus'],
                     stimData[, "UASEvents"],
                     stimData[, which(colnames(stimData)=="UASLAeq"):
                                which(colnames(stimData)=="SNRlevel")],
                     stimData[, which(colnames(stimData)=="IntermitRatioC2MaxLR"):
                                which(colnames(stimData)=="IntermitRatioC5MaxLR")],
                     stimData[, which(colnames(stimData)=="UASLAEMaxLR"):
                                which(colnames(stimData)=="UASEPNLMaxLR")],
                     stimData[, which(colnames(stimData)=="UASLoudECMAPowAvgBin"):
                                which(colnames(stimData)=="UASLoudISO3PowAvgBin")],
                     stimData[, which(colnames(stimData)=="UASTonalECMAAvgMaxLR"):
                                which(colnames(stimData)=="UASSharpvBISO105ExBin")],
                     stimData[, which(colnames(stimData)=="UASImpulsSHMPowAvgMaxLR"):
                                which(colnames(stimData)=="UASPsychAnnoyBoucher")],
                     stimData[, which(colnames(stimData)=="LAeqLAF90diff"):
                                which(colnames(stimData)=="dPsychAnnoyBoucher")],
                     stimData[, which(colnames(stimData)=="ValenceMedian"):
                                which(colnames(stimData)=="dHighAnnoyProp")],
                     stimData[, which(colnames(stimData)=="HighAnnoyPc"):
                                which(colnames(stimData)=="dHighAnnoyPc")])

# remove duplicated variables
stimDataNum <- subset(stimDataNum, select = -c(UASLAeq))

colnames(stimDataNum)[1] <- "Stimulus"
colnames(stimDataNum)[2] <- "UASEvents"

# make the discrete ordinal outcome variables factors
stimDataNum[['UASEvents']] <- factor(stimDataNum[['UASEvents']], levels=c(0, 1, 3, 5, 9), order=TRUE)
stimDataNum[['ValenceMedian']] <- factor(stimDataNum[['ValenceMedian']], levels=c(1, 1.5, 2, 2.5, 3, 3.5, 4, 4.5, 5), order=TRUE)
stimDataNum[['ArousalMedian']] <- factor(stimDataNum[['ArousalMedian']], levels=c(1, 1.5, 2, 2.5, 3, 3.5, 4, 4.5, 5), order=TRUE)
stimDataNum[['AnnoyMedian']] <- factor(stimDataNum[['AnnoyMedian']], levels=c(0, 0.5, 1, 1.5, 2, 2.5, 3, 3.5, 4, 4.5, 5,
                                                                                5.5, 6, 6.5, 7, 7.5, 8, 8.5, 9, 9.5, 10), order=TRUE)
stimDataNum[['dValenceMedian']] <- factor(stimDataNum[['dValenceMedian']], levels=c(-4, -3.5, -3, -2.5, -2, -1.5, -1, -0.5, 0,
                                                                                      0.5, 1, 1.5, 2, 2.5, 3, 3.5,  4), order=TRUE)
stimDataNum[['dArousalMedian']] <- factor(stimDataNum[['dArousalMedian']], levels=c(-4, -3.5, -3, -2.5, -2, -1.5, -1, -0.5, 0,
                                                                                      0.5, 1, 1.5, 2, 2.5, 3, 3.5,  4), order=TRUE)
stimDataNum[['dAnnoyMedian']] <- factor(stimDataNum[['dAnnoyMedian']], levels=c(-10, -9.5, -9, -8.5, -8, -7.5, -7, -6.5, -6, -5.5, -5,
                                                                                  -4.5, -4, -3.5, -3, -2.5, -2, -1.5, -1, -0.5,
                                                                                 0, 0.5, 1, 1.5, 2, 2.5, 3, 3.5, 4, 4.5, 5,
                                                                                 5.5, 6, 6.5, 7, 7.5, 8, 8.5, 9, 9.5, 10), order=TRUE)

# omit ambient-only stimuli
stimDataNum <- stimDataNum |> dplyr::filter(UASEvents != 0)


stimDataNum$SNRlevel <- as.numeric(stimDataNum$SNRlevel)

Random forest functions

Write a function to train a conditional-inference random forest (crf) model on input data according to input formula, iterate over input random seeds, average error and variable importance metrics, and output metrics with plotted

Averaging over multiple random seeds

multi_crfReg <- function(dataIn, iVars, dVar, seeds, ntree, mtry, permImpCondThres=0.95, minsplit=20, minbucket=7, nperm=1){
  # initialise variables
  crfOOBErrAll <- 0
  crfOOBRMSE <- 0
  crfOOBMAE <- 0
  crfOOBErrR2 <- 0
  crfMarPermImpVals <- 0
  crfConPermImpVals <- 0
  crfMarPermImpValsPerTree <- data.frame()
  crfConPermImpValsPerTree <- data.frame()
  
  for (iters in 1:length(seeds)){
    
    # formula for regression
    formVars <- reformulate(iVars, dVar)
    
    # set random seed
    set.seed(seeds[iters])
    # train crf model
    crfModel <- party::cforest(formVars, data=dataIn,
                               controls=party::cforest_unbiased(ntree=ntree,
                                                                mtry=mtry,
                                                                minsplit=minsplit,
                                                                minbucket=minbucket))
    
    # get OOB predictions
    crfModelOOB <- predict(crfModel, OOB=TRUE, type='response')
    
    # get OOB error
    crfModelOOBErr <- as.numeric(as.matrix(as.numeric(as.matrix(crfModelOOB))
                                           - as.numeric(as.matrix(crfModel@data@env$response[[names(crfModel@data@env$response)]]))))

    # OOB RMSE, MAE and Rsquared
    crfOOBRMSE <- crfOOBRMSE + sqrt(mean(crfModelOOBErr^2))
    crfOOBMAE <- crfOOBMAE + mean(abs(crfModelOOBErr))
    crfOOBErrR2 <- crfOOBErrR2 + cor(as.numeric(as.matrix(crfModelOOB)),
                                     as.numeric(as.matrix(crfModel@data@env$response[[names(crfModel@data@env$response)]])))^2

    # set random seed
    set.seed(seeds[iters])

    # set random seed
    set.seed(seeds[iters])
    # conditional variable permutation importance
    crfConPermImp <- permimp::permimp(crfModel, nperm=nperm, conditional=TRUE, threshold=permImpCondThres, progressBar=FALSE)
    
    crfConPermImpVals <- crfConPermImpVals + crfConPermImp$values
    crfConPermImpValsPerTree <- rbind(crfConPermImpValsPerTree, crfConPermImp$perTree)
  }
  
  # average metrics
  crfOOBErrAll <- crfOOBErrAll/length(seeds)
  crfOOBRMSE <- crfOOBRMSE/length(seeds)
  crfOOBMAE <- crfOOBMAE/length(seeds)
  crfOOBErrR2 <- crfOOBErrR2/length(seeds)
  crfConPermImpVals <- data.frame(CondPermImp=sort(crfConPermImpVals/length(seeds), decreasing=TRUE))
  crfConPermImpValsQtl <- data.frame(apply(crfConPermImpValsPerTree, 2, quantile, probs=c(0.25, 0.50, 0.75)))
  
  resultsOut <- list('OOB_RMSE'=crfOOBRMSE, 'OOB_MAE'=crfOOBMAE, 'Rsquared'=crfOOBErrR2, 'conditional_permimp'=crfConPermImpVals,                      'conditional_permimp_perTree'=crfConPermImpValsPerTree, 'conditional_permimp_qtl'=crfConPermImpValsQtl)
  return(resultsOut)
}

Comparing rankings from two seeds

crfReg <- function(dataIn, iVars, dVar, seeds, ntree, mtry, permImpCondThres=0.95, minsplit=20, minbucket=7, nperm=1){
  # initialise variables
  crfOOBErrAll <- 0
  crfOOBRMSE <- 0
  crfOOBMAE <- 0
  crfOOBErrR2 <- 0
  crfMarPermImpVals <- 0
  crfConPermImpVals <- 0
  crfMarPermImpValsPerTree <- data.frame()
  crfConPermImpValsPerTree <- data.frame()

  # formula for regression
  formVars <- reformulate(iVars, dVar)
  
  for (iters in 1:length(seeds)){
  
    # set random seed
    set.seed(seeds[iters])
    # train crf model
    crfModel <- party::cforest(formVars, data=dataIn,
                               controls=party::cforest_unbiased(ntree=ntree,
                                                                mtry=mtry,
                                                                minsplit=minsplit,
                                                                minbucket=minbucket))
    
    # conditional variable permutation importance
    crfConPermImp <- permimp::permimp(crfModel, nperm=nperm, conditional=TRUE, threshold=permImpCondThres, progressBar=FALSE)
    
    crfConPermImpVals <- crfConPermImp$values
    
    if (iters == 1){
      crfConPermImpVals1 <- data.frame(CondPermImp=sort(crfConPermImpVals, decreasing=TRUE))
      crfConPermImpValsPerTree1 <- crfConPermImp$perTree
      crfConPermImpValsQtl1 <- data.frame(apply(crfConPermImpValsPerTree1, 2, quantile, probs=c(0.25, 0.50, 0.75)))
      
      # get OOB predictions
      crfModelOOB <- predict(crfModel, OOB=TRUE, type='response')
      
      # get OOB error
      crfModelOOBErr <- as.numeric(as.matrix(as.numeric(as.matrix(crfModelOOB))
                                              - as.numeric(as.matrix(crfModel@data@env$response[[names(crfModel@data@env$response)]]))))
      
      # OOB RMSE, error quartiles and Rsquared
      crfOOBRMSE <- sqrt(mean(crfModelOOBErr^2))
      crfOOBMAE <- crfOOBMAE + mean(abs(crfModelOOBErr))
      crfOOBErrR2 <- cor(as.numeric(as.matrix(crfModelOOB)),
                                    as.numeric(as.matrix(crfModel@data@env$response[[names(crfModel@data@env$response)]])))^2

      }
    
    else{
      crfConPermImpValsN <- data.frame(CondPermImp=sort(crfConPermImpVals, decreasing=TRUE))
      
      nVarImpChecks <- min(max(sum(crfConPermImpVals1 >= crfConPermImpVals1$CondPermImp[1]*0.1),
                               sum(crfConPermImpValsN >= crfConPermImpValsN$CondPermImp[1]*0.1)), 4)  # number of variable importance values with a value at least 10% of the highest importance
      if (sum(rownames(crfConPermImpVals1)[1:nVarImpChecks] != rownames(crfConPermImpValsN)[1:nVarImpChecks]) > 0){
        warning("Permutation importance rank order within 10% of max differs between seeds: increase number of trees ('ntree') or number of permutations ('nperm'), or subsample of features ('mtry')")
      }
      else{
        resultsOut <- list('OOB_errors'=crfModelOOBErr, 'OOB_RMSE'=crfOOBRMSE, 'OOB_MAE'=crfOOBMAE, 'Rsquared'=crfOOBErrR2, 'conditional_permimp'=crfConPermImpVals1, 'conditional_permimp_perTree'=crfConPermImpValsPerTree1, 'conditional_permpimp_qtl'=crfConPermImpValsQtl1)
        return(resultsOut)
      }
      
    }
    
  }

}

Hyperparameter tuning

length(iVars)

[1] 234

Parts A & B analysis

Set global parameters

permImpCondThres <- 0.95
minsplit <- 20
minbucket <- 7
ntrees <- c(251, 501, 1001, 1501, 2501, 4001, 5501)

eventVar <- "UASEvents"
ambVar <- "AmbientLAeq"

Mean change in annoyance

Initialise results output variables

resdAnnoyMnFitAB <- data.frame(RMSE = numeric(),
                               MAE = numeric(),
                               Rsquared = numeric())
resdAnnoyMnPermImpAB <- list()

All variables (absolute and difference)

Set variables

iVars <- names(stimDataNum)[which(names(stimDataNum) == 'UASEvents'):which(names(stimDataNum) == 'dPsychAnnoyBoucher')]

dVar <- "dAnnoyMean"

seeds <- c(14569, 98651, 54654498, 454948, 41321)

Hyperparameter tuning

p <- mtryTune(dataIn=stimDataNum, iVars=iVars, dVar=dVar, seeds=seeds[1:2],
              ntrees=ntrees, minsplit=minsplit, minbucket=minbucket)

p

if (saveplots){
  ggsave(filename="PtsABdAnnoyMnAllVarsHyperTune.svg", width=12, height=4, path=file.path(outFigPath, "svg"))
  unlink("PtsABdAnnoyMnAllVarsHyperTune.svg")

  ggsave(filename="PtsABdAnnoyMnAllVarsHyperTune.pdf", width=12, height=4, path=file.path(outFigPath, "pdf"))
  unlink("PtsABdAnnoyMnAllVarsHyperTune.pdf")
}

Selected hyperparameters

ntree <- 2501
mtry <- as.integer(length(iVars)/3.5)

Run model

Train preliminary model

nperm <- 5

resultsOutAbsDiffs <- crfReg(dataIn=stimDataNum, iVars=iVars, dVar=dVar, seeds=seeds[1:2], ntree=ntree, mtry=mtry, permImpCondThres=permImpCondThres, nperm=nperm, minsplit=minsplit, minbucket=minbucket)

# print model prediction results
resultsOutAbsDiffs$OOB_RMSE

[1] 0.5232516

resultsOutAbsDiffs$OOB_MAE

[1] 0.4044405

resultsOutAbsDiffs$Rsquared

[1] 0.8933325

Train multiple seeds model

resultsOutAbsDiffs <- multi_crfReg(dataIn=stimDataNum, iVars=iVars, dVar=dVar, seeds=seeds, ntree=ntree, mtry=mtry, permImpCondThres=permImpCondThres, nperm=nperm, minsplit=minsplit, minbucket=minbucket)

# print model prediction results
resultsOutAbsDiffs$OOB_RMSE

[1] 0.5251036

resultsOutAbsDiffs$OOB_MAE

[1] 0.4057256

resultsOutAbsDiffs$Rsquared

[1] 0.8925424

# store results
resdAnnoyMnFitAB['All vars', 'RMSE'] <- resultsOutAbsDiffs$OOB_RMSE
resdAnnoyMnFitAB['All vars', 'MAE'] <- resultsOutAbsDiffs$OOB_MAE
resdAnnoyMnFitAB['All vars', 'Rsquared'] <- resultsOutAbsDiffs$Rsquared
resdAnnoyMnPermImpAB$AllVars <- resultsOutAbsDiffs$conditional_permimp

Plot results

par(mai=c(0,3,0,0))

# plot conditional importance
resultsOutAbsDiffs.conimp <- arrange(resultsOutAbsDiffs$conditional_permimp, desc(row_number()))

pBar <- ggplot(resultsOutAbsDiffs.conimp) + geom_col(aes(x=factor(rownames(resultsOutAbsDiffs.conimp), levels=rownames(resultsOutAbsDiffs.conimp)), y=CondPermImp), fill=mycolours[9], width=0.5) + labs(x="Variable", y="Conditional variable permutation importance (mean change in annoyance)") + theme(text = element_text(family = "serif"), panel.grid=element_line(color = rgb(235, 235, 235, 100, maxColorValue = 255), linewidth = 0.25, linetype = 2)) + coord_flip()
pBar

if (saveplots){
  ggsave(filename="PtsABdAnnoyMnAllVarsConPermimp.svg", width=8, height=30, path=file.path(outFigPath, "svg"))
  unlink("PtsABdAnnoyMnAllVarsConPermimp.svg")

  ggsave(filename="PtsABdAnnoyMnAllVarsConPermimp.pdf", width=8, height=30, path=file.path(outFigPath, "pdf"))
  unlink("PtsABdAnnoyMnAllVarsConPermimp.pdf")
}

# Plot only positive values

resultsOutAbsDiffs.conimpPtv <- resultsOutAbsDiffs.conimp |>
                                          rownames_to_column('Metric') |>
                                                filter_if(is.numeric, all_vars(. > 0)) |>
                                                      column_to_rownames('Metric')

pBar <- ggplot(resultsOutAbsDiffs.conimpPtv) + geom_col(aes(x=factor(rownames(resultsOutAbsDiffs.conimpPtv), levels=rownames(resultsOutAbsDiffs.conimpPtv)), y=CondPermImp), fill=mycolours[9], width=0.5) + labs(x="Variable", y="Conditional variable permutation importance (mean change in annoyance)") + theme(text = element_text(family = "serif"), panel.grid=element_line(color = rgb(235, 235, 235, 100, maxColorValue = 255), linewidth = 0.25, linetype = 2)) + coord_flip()
pBar

if (saveplots){
  ggsave(filename="PtsABdAnnoyMnAllVarsConPermimpPtv.svg", width=8, height=18, path=file.path(outFigPath, "svg"))
  unlink("PtsABdAnnoyMnAllVarsConPermimp.svg")
  
  ggsave(filename="PtsABdAnnoyMnAllVarsConPermimpPtv.pdf", width=8, height=18, path=file.path(outFigPath, "pdf"))
  unlink("PtsABdAnnoyMnAllVarsConPermimp.pdf")
}

# Plot only values within 1% of the maximum

resultsOutAbsDiffs.conimp1pc <- resultsOutAbsDiffs.conimp |>
                                          rownames_to_column('Metric') |>
                                                filter_if(is.numeric, all_vars(. > max(resultsOutAbsDiffs.conimp)/100)) |>
                                                      column_to_rownames('Metric')

pBar <- ggplot(resultsOutAbsDiffs.conimp1pc) + geom_col(aes(x=factor(rownames(resultsOutAbsDiffs.conimp1pc), levels=rownames(resultsOutAbsDiffs.conimp1pc)), y=CondPermImp), fill=mycolours[9], width=0.5) + labs(x="Variable", y="Conditional variable permutation importance (mean change in annoyance)") + theme(text = element_text(family = "serif"), panel.grid=element_line(color = rgb(235, 235, 235, 100, maxColorValue = 255), linewidth = 0.25, linetype = 2)) + coord_flip()
pBar

if (saveplots){
  ggsave(filename="PtsABdAnnoyMnAllVarsConPermimp1pc.svg", width=8, height=6, path=file.path(outFigPath, "svg"))
  unlink("PtsABdAnnoyMnAllVarsConPermimp1pc.svg")
  
  ggsave(filename="PtsABdAnnoyMnAllVarsConPermimp1pc.pdf", width=8, height=6, path=file.path(outFigPath, "pdf"))
  unlink("PtsABdAnnoyMnAllVarsConPermimp1pc.pdf")
}

Absolute variables

Set variables

iVars <- names(stimDataNum)[which(names(stimDataNum) == 'UASEvents'):which(names(stimDataNum) == 'UASPsychAnnoyBoucher')]
iVars <- iVars[! iVars %in% c('SNRlevel', 'IntermitRatioC2MaxLR', 'IntermitRatioC3MaxLR', 'IntermitRatioC5MaxLR')]

dVar <- "dAnnoyMean"

seeds <- c(578312, 544, 84894, 54654, 153157)

Hyperparameter tuning

p <- mtryTune(dataIn=stimDataNum, iVars=iVars, dVar=dVar, seeds=seeds[1:2],
              ntrees=ntrees, minsplit=minsplit, minbucket=minbucket)

p

if (saveplots){
  ggsave(filename="PtsABdAnnoyMnAbsVarsHyperTune.svg", width=12, height=4, path=file.path(outFigPath, "svg"))
  unlink("PtsABdAnnoyMnAbsVarsHyperTune.svg")

  ggsave(filename="PtsABdAnnoyMnAbsVarsHyperTune.pdf", width=12, height=4, path=file.path(outFigPath, "pdf"))
  unlink("PtsABdAnnoyMnAbsVarsHyperTune.pdf")
}

Selected hyperparameters

ntree <- 2501
mtry <- as.integer(length(iVars)/1.75)

Run model

Train preliminary model

nperm <- 5

resultsOutAbs <- crfReg(dataIn=stimDataNum, iVars=iVars, dVar=dVar, seeds=seeds[1:2], ntree=ntree, mtry=mtry, permImpCondThres=permImpCondThres, nperm=nperm, minsplit=minsplit, minbucket=minbucket)

# print model prediction results
resultsOutAbs$OOB_RMSE

[1] 0.7179463

resultsOutAbs$OOB_MAE

[1] 0.5697517

resultsOutAbs$Rsquared

[1] 0.81913

Train multiple seeds model

resultsOutAbs <- multi_crfReg(dataIn=stimDataNum, iVars=iVars, dVar=dVar, seeds=seeds, ntree=ntree, mtry=mtry, permImpCondThres=permImpCondThres, nperm=nperm, minsplit=minsplit, minbucket=minbucket)

# print model prediction results
resultsOutAbs$OOB_RMSE

[1] 0.7262005

resultsOutAbs$OOB_MAE

[1] 0.5771998

resultsOutAbs$Rsquared

[1] 0.8140447

# store results
resdAnnoyMnFitAB['Abs vars', 'RMSE'] <- resultsOutAbs$OOB_RMSE
resdAnnoyMnFitAB['Abs vars', 'MAE'] <- resultsOutAbs$OOB_MAE
resdAnnoyMnFitAB['Abs vars', 'Rsquared'] <- resultsOutAbs$Rsquared
resdAnnoyMnPermImpAB$AbsVars <- resultsOutAbs$conditional_permimp

Plot results

par(mai=c(0,3,0,0))

# plot conditional importance
resultsOutAbs.conimp <- arrange(resultsOutAbs$conditional_permimp, desc(row_number()))

pBar <- ggplot(resultsOutAbs.conimp) + geom_col(aes(x=factor(rownames(resultsOutAbs.conimp), levels=rownames(resultsOutAbs.conimp)), y=CondPermImp), fill=mycolours[1], width=0.5) + labs(x="Variable", y="Conditional variable permutation importance (mean change in annoyance)") + theme(text = element_text(family = "serif"), panel.grid=element_line(color = rgb(235, 235, 235, 100, maxColorValue = 255), linewidth = 0.25, linetype = 2)) +
  coord_flip()
pBar

if (saveplots){
  ggsave(filename="PtsABdAnnoyMnAbsVarsConPermimp.svg", width=8, height=13, path=file.path(outFigPath, "svg"))
  unlink("PtsABdAnnoyMnAbsVarsConPermimp.svg")
  
  ggsave(filename="PtsABdAnnoyMnAbsVarsConPermimp.pdf", width=8, height=13, path=file.path(outFigPath, "pdf"))
  unlink("PtsABdAnnoyMnAbsVarsConPermimp.pdf")
}

# Plot only positive values

resultsOutAbs.conimpPtv <- resultsOutAbs.conimp |>
                                          rownames_to_column('Metric') |>
                                                filter_if(is.numeric, all_vars(. > 0)) |>
                                                      column_to_rownames('Metric')

pBar <- ggplot(resultsOutAbs.conimpPtv) + geom_col(aes(x=factor(rownames(resultsOutAbs.conimpPtv), levels=rownames(resultsOutAbs.conimpPtv)), y=CondPermImp), fill=mycolours[1], width=0.5) + labs(x="Variable", y="Conditional variable permutation importance (mean change in annoyance)") + theme(text = element_text(family = "serif"), panel.grid=element_line(color = rgb(235, 235, 235, 100, maxColorValue = 255), linewidth = 0.25, linetype = 2)) + coord_flip()
pBar

if (saveplots){
  ggsave(filename="PtsABdAnnoyMnAbsVarsConPermimpPtv.svg", width=8, height=10, path=file.path(outFigPath, "svg"))
  unlink("PtsABdAnnoyMnAbsVarsConPermimpPtv.svg")
  
  ggsave(filename="PtsABdAnnoyMnAbsVarsConPermimpPtv.pdf", width=8, height=10, path=file.path(outFigPath, "pdf"))
  unlink("PtsABdAnnoyMnAbsVarsConPermimpPtv.pdf")
}

# Plot only values within 1% of the maximum

resultsOutAbs.conimp1pc <- resultsOutAbs.conimp |>
                                          rownames_to_column('Metric') |>
                                                filter_if(is.numeric, all_vars(. > max(resultsOutAbs.conimp)/100)) |>
                                                      column_to_rownames('Metric')

pBar <- ggplot(resultsOutAbs.conimp1pc) + geom_col(aes(x=factor(rownames(resultsOutAbs.conimp1pc), levels=rownames(resultsOutAbs.conimp1pc)), y=CondPermImp), fill=mycolours[1], width=0.5) + labs(x="Variable", y="Conditional variable permutation importance (mean change in annoyance)") + theme(text = element_text(family = "serif"), panel.grid=element_line(color = rgb(235, 235, 235, 100, maxColorValue = 255), linewidth = 0.25, linetype = 2)) + coord_flip()
pBar

if (saveplots){
  ggsave(filename="PtsABdAnnoyMnAbsVarsConPermimp1pc.svg", width=8, height=3, path=file.path(outFigPath, "svg"))
  unlink("PtsABdAnnoyMnAbsVarsConPermimp1pc.svg")
  
  ggsave(filename="PtsABdAnnoyMnAbsVarsConPermimp1pc.pdf", width=8, height=3, path=file.path(outFigPath, "pdf"))
  unlink("PtsABdAnnoyMnAbsVarsConPermimp1pc.pdf")
}

Selected metric

absVar <- "UASLAEMaxLR"

SQM analysis

Individual SQMs

Sharpness

Set variables

iVars <- c(absVar, eventVar, ambVar, "UASSharpAurISO3PowAvgBin", "UASSharpAurISO305ExBin", "UASSharpAurSHMPowAvgBin", "UASSharpAurSHM05ExBin", "UASSharpAurISO1PowAvgBin", "UASSharpAurISO105ExBin", "UASSharpvBISO1PowAvgBin", "UASSharpvBISO105ExBin", "UASSharpDINPowAvgBin", "UASSharpDIN05ExBin", "UASSharpAurISO1MedBin",
         "UASTonShpAurSHMPowAvgBin", "UASTonShpAurSHM05ExBin")
dVar <- "dAnnoyMean"

seeds <- c(7041, 905, 4984651, 6513213, 120651)

Hyperparameter tuning

p <- mtryTune(dataIn=stimDataNum, iVars=iVars, dVar=dVar, seeds=seeds[1:2],
             ntrees=ntrees, minsplit=minsplit, minbucket=minbucket)

p

Selected hyperparameters

ntree <- 251
mtry <- as.integer(length(iVars)/1.5)

Run model

Train preliminary model

nperm <- 5

resultsOutSharp <- crfReg(dataIn=stimDataNum, iVars=iVars, dVar=dVar, seeds=seeds[1:2], ntree=ntree, mtry=mtry, permImpCondThres=permImpCondThres, nperm=nperm, minsplit=minsplit, minbucket=minbucket)

# print model prediction results
resultsOutSharp$OOB_RMSE

[1] 0.6520071

resultsOutSharp$OOB_MAE

[1] 0.5246782

resultsOutSharp$Rsquared

[1] 0.8553724

Train multiple seeds model

resultsOutSharp <- multi_crfReg(dataIn=stimDataNum, iVars=iVars, dVar=dVar, seeds=seeds, ntree=ntree, mtry=mtry, permImpCondThres=permImpCondThres, nperm=nperm, minsplit=minsplit, minbucket=minbucket)

# print model prediction results
resultsOutSharp$OOB_RMSE

[1] 0.6548829

resultsOutSharp$OOB_MAE

[1] 0.5229297

resultsOutSharp$Rsquared

[1] 0.8539235

# store results
resdAnnoyMnFitAB['Abs sharp', 'RMSE'] <- resultsOutSharp$OOB_RMSE
resdAnnoyMnFitAB['Abs sharp', 'MAE'] <- resultsOutSharp$OOB_MAE
resdAnnoyMnFitAB['Abs sharp', 'Rsquared'] <- resultsOutSharp$Rsquared
resdAnnoyMnPermImpAB$AbsSharp <- resultsOutSharp$conditional_permimp

Plot results

par(mai=c(0,3,0,0))

# plot conditional importance
resultsOutSharp.conimp <- arrange(resultsOutSharp$conditional_permimp, desc(row_number()))

pBar <- ggplot(resultsOutSharp.conimp) + geom_col(aes(x=factor(rownames(resultsOutSharp.conimp), levels=rownames(resultsOutSharp.conimp)), y=CondPermImp), fill=mycolours[2], width=0.5) + labs(x="Variable", y="Conditional variable permutation importance (mean change in annoyance)") + ggtitle("Sharpness") + theme(text = element_text(family = "serif"), panel.grid=element_line(color = rgb(235, 235, 235, 100, maxColorValue = 255), linewidth = 0.25, linetype = 2)) + coord_flip()
pBar

if (saveplots){
  ggsave(filename="PtsABdAnnoyMnSharpConPermimp.svg", width=8, height=4.9, path=file.path(outFigPath, "svg"))
  unlink("PtsABdAnnoyMnSharpConPermimp.svg")
  
  ggsave(filename="PtsABdAnnoyMnSharpConPermimp.pdf", width=8, height=4.9, path=file.path(outFigPath, "pdf"))
  unlink("PtsABdAnnoyMnSharpConPermimp.pdf")
}

Selected metric

sharpVar <- "UASSharpAurISO305ExBin"

Tonal loudness and tonality

Set variables

iVars <- c(absVar, eventVar, ambVar, "UASTonalECMAAvgMaxLR", "UASTonalSHMInt05ExMaxLR", "UASTonalSHMIntAvgMaxLR", "UASTonalECMA05ExMaxLR", "UASTonalAwSHMAvgMaxLR", "UASTonalAwSHM05ExMaxLR",   "UASTonalAwSHMIntAvgMaxLR",     "UASTonalAwSHMInt05ExMaxLR", "UASTonLdECMAPowAvgBin", "UASTonLdECMA05ExBin", "UASTonalAurAvgMaxLR", "UASTonalAur05ExMaxLR", "UASTonalAur10ExMaxLR",
         "UASTonShpAurSHMPowAvgBin", "UASTonShpAurSHM05ExBin")
dVar <- "dAnnoyMean"

seeds <- c(540, 104798, 456464, 87331, 94564)

Hyperparameter tuning

p <- mtryTune(dataIn=stimDataNum, iVars=iVars, dVar=dVar, seeds=seeds[1:2],
             ntrees=ntrees, minsplit=minsplit, minbucket=minbucket)

p

Selected hyperparameters

ntree <- 251
mtry <- as.integer(length(iVars)/1.5)

Run model

Train preliminary model

# Tonality with tonal loudness

nperm <- 5

resultsOutTonal1 <- crfReg(dataIn=stimDataNum, iVars=iVars, dVar=dVar, seeds=seeds[1:2], ntree=ntree, mtry=mtry, permImpCondThres=permImpCondThres, nperm=nperm, minsplit=minsplit, minbucket=minbucket)

# print model prediction results
resultsOutTonal1$OOB_RMSE

[1] 0.6817226

resultsOutTonal1$OOB_MAE

[1] 0.5354283

resultsOutTonal1$Rsquared

[1] 0.8328264

Train multiple seeds model

# Tonality with tonal loudness

resultsOutTonal1 <- multi_crfReg(dataIn=stimDataNum, iVars=iVars, dVar=dVar, seeds=seeds, ntree=ntree, mtry=mtry, permImpCondThres=permImpCondThres, nperm=nperm, minsplit=minsplit, minbucket=minbucket)

# print model prediction results
resultsOutTonal1$OOB_RMSE

[1] 0.6726346

resultsOutTonal1$OOB_MAE

[1] 0.5320356

resultsOutTonal1$Rsquared

[1] 0.8402319

# store results
resdAnnoyMnFitAB['Abs tonal inc loud', 'RMSE'] <- resultsOutTonal1$OOB_RMSE
resdAnnoyMnFitAB['Abs tonal inc loud', 'MAE'] <- resultsOutTonal1$OOB_MAE
resdAnnoyMnFitAB['Abs tonal inc loud', 'Rsquared'] <- resultsOutTonal1$Rsquared
resdAnnoyMnPermImpAB$AbsTonal1 <- resultsOutTonal1$conditional_permimp

Plot results

par(mai=c(0,3,0,0))

# plot conditional importance
resultsOutTonal1.conimp <- arrange(resultsOutTonal1$conditional_permimp, desc(row_number()))

pBar <- ggplot(resultsOutTonal1.conimp) + geom_col(aes(x=factor(rownames(resultsOutTonal1.conimp), levels=rownames(resultsOutTonal1.conimp)), y=CondPermImp), fill=mycolours[3], width=0.5) + labs(x="Variable", y="Conditional variable permutation importance (mean change in annoyance)") + ggtitle("Tonality inc. tonal loudness") + theme(text = element_text(family = "serif"), panel.grid=element_line(color = rgb(235, 235, 235, 100, maxColorValue = 255), linewidth = 0.25, linetype = 2)) + coord_flip(ylim=c(0, 1.4))
pBar

if (saveplots){
  ggsave(filename="PtsABdAnnoyMnTonalLdConPermimp.svg", width=8, height=5, path=file.path(outFigPath, "svg"))
  unlink("PtsABdAnnoyMnTonalLdConPermimp.svg")
  
  ggsave(filename="PtsABdAnnoyMnTonalLdConPermimp.pdf", width=8, height=5, path=file.path(outFigPath, "pdf"))
  unlink("PtsABdAnnoyMnTonalLdConPermimp.pdf")
}

Selected metric

tonLdVar <- "UASTonLdECMAPowAvgBin"

Tonality without tonal loudness or tonal sharpness

Set variables

iVars <- c(absVar, eventVar, ambVar, "UASTonalECMAAvgMaxLR", "UASTonalSHMInt05ExMaxLR", "UASTonalSHMIntAvgMaxLR", "UASTonalECMA05ExMaxLR", "UASTonalAwSHMAvgMaxLR", "UASTonalAwSHM05ExMaxLR",   "UASTonalAwSHMIntAvgMaxLR", "UASTonalAwSHMInt05ExMaxLR",    "UASTonalAurAvgMaxLR", "UASTonalAur05ExMaxLR", "UASTonalAur10ExMaxLR")
dVar <- "dAnnoyMean"

seeds <- c(156089, 5860, 10528, 89541, 4685146)

Hyperparameter tuning

p <- mtryTune(dataIn=stimDataNum, iVars=iVars, dVar=dVar, seeds=seeds[1:2],
             ntrees=ntrees, minsplit=minsplit, minbucket=minbucket)

p

Selected hyperparameters

ntree <- 501
mtry <- as.integer(length(iVars)/1.25)

Run model

Train preliminary model

# Tonality

nperm <- 5

resultsOutTonal2 <- crfReg(dataIn=stimDataNum, iVars=iVars, dVar=dVar, seeds=seeds[1:2], ntree=ntree, mtry=mtry, permImpCondThres=permImpCondThres, nperm=nperm, minsplit=minsplit, minbucket=minbucket)

# print model prediction results
resultsOutTonal2$OOB_RMSE

[1] 0.6735604

resultsOutTonal2$OOB_MAE

[1] 0.5215824

resultsOutTonal2$Rsquared

[1] 0.835001

Train multiple seeds model

# Tonality

resultsOutTonal2 <- multi_crfReg(dataIn=stimDataNum, iVars=iVars, dVar=dVar, seeds=seeds, ntree=ntree, mtry=mtry, permImpCondThres=permImpCondThres, nperm=nperm, minsplit=minsplit, minbucket=minbucket)

# print model prediction results
resultsOutTonal2$OOB_RMSE

[1] 0.6770415

resultsOutTonal2$OOB_MAE

[1] 0.526503

resultsOutTonal2$Rsquared

[1] 0.8322548

# store results
resdAnnoyMnFitAB['Abs tonal no loud', 'RMSE'] <- resultsOutTonal2$OOB_RMSE
resdAnnoyMnFitAB['Abs tonal no loud', 'MAE'] <- resultsOutTonal2$OOB_MAE
resdAnnoyMnFitAB['Abs tonal no loud', 'Rsquared'] <- resultsOutTonal2$Rsquared
resdAnnoyMnPermImpAB$AbsTonal2 <- resultsOutTonal2$conditional_permimp

Plot results

par(mai=c(0,3,0,0))

# plot conditional importance
resultsOutTonal2.conimp <- arrange(resultsOutTonal2$conditional_permimp, desc(row_number()))

pBar <- ggplot(resultsOutTonal2.conimp) + geom_col(aes(x=factor(rownames(resultsOutTonal2.conimp), levels=rownames(resultsOutTonal2.conimp)), y=CondPermImp), fill=mycolours[3], width=0.5) + labs(x="Variable", y="Conditional variable permutation importance (mean change in annoyance)") + ggtitle("Tonality w/o tonal loudness") + theme(text = element_text(family = "serif"), panel.grid=element_line(color = rgb(235, 235, 235, 100, maxColorValue = 255), linewidth = 0.25, linetype = 2)) + coord_flip(ylim=c(0, 1.4))
pBar

if (saveplots){
  ggsave(filename="PtsABdAnnoyMnTonalConPermimp.svg", width=8, height=4.4, path=file.path(outFigPath, "svg"))
  unlink("PtsABdAnnoyMnTonalConPermimp.svg")
  
  ggsave(filename="PtsABdAnnoyMnTonalConPermimp.pdf", width=8, height=4.4, path=file.path(outFigPath, "pdf"))
  unlink("PtsABdAnnoyMnTonalConPermimp.pdf")
}

Selected metric

tonalVar <- "UASTonalAwSHMInt05ExMaxLR"

Fluctuation strength

Set variables

# Fluctuation strength
iVars <- c(absVar, eventVar, ambVar, "UASFluctOldSHM10ExBin", "UASFluctOldSHM05ExBin", "UASFluctECMA10ExBin", "UASFluctECMA05ExBin", "UASFluctFZ10ExMaxLR", "UASFluctFZ05ExMaxLR", "UASFluctOV10ExMaxLR", "UASFluctOV05ExMaxLR")
dVar <- "dAnnoyMean"

seeds <- c(25107, 546098, 195, 5937, 102658)

Hyperparameter tuning

p <- mtryTune(dataIn=stimDataNum, iVars=iVars, dVar=dVar, seeds=seeds[1:2],
             ntrees=ntrees, minsplit=minsplit, minbucket=minbucket)

p

Selected hyperparameters

ntree <- 5501
mtry <- as.integer(length(iVars)/1.5)

Run model

Train preliminary model

nperm <- 5

resultsOutFluct <- crfReg(dataIn=stimDataNum, iVars=iVars, dVar=dVar, seeds=seeds[1:2], ntree=ntree, mtry=mtry, permImpCondThres=permImpCondThres, nperm=nperm, minsplit=minsplit, minbucket=minbucket)

# print model prediction results
resultsOutFluct$OOB_RMSE

[1] 0.6411068

resultsOutFluct$OOB_MAE

[1] 0.5139778

resultsOutFluct$Rsquared

[1] 0.8590107

Train multiple seeds model

resultsOutFluct <- multi_crfReg(dataIn=stimDataNum, iVars=iVars, dVar=dVar, seeds=seeds, ntree=ntree, mtry=mtry, permImpCondThres=permImpCondThres, nperm=nperm, minsplit=minsplit, minbucket=minbucket)

# print model prediction results
resultsOutFluct$OOB_RMSE

[1] 0.6442208

resultsOutFluct$OOB_MAE

[1] 0.5160342

resultsOutFluct$Rsquared

[1] 0.8570714

# store results
resdAnnoyMnFitAB['Abs fluct', 'RMSE'] <- resultsOutFluct$OOB_RMSE
resdAnnoyMnFitAB['Abs fluct', 'MAE'] <- resultsOutFluct$OOB_MAE
resdAnnoyMnFitAB['Abs fluct', 'Rsquared'] <- resultsOutFluct$Rsquared
resdAnnoyMnPermImpAB$AbsFluct <- resultsOutFluct$conditional_permimp

Plot results

par(mai=c(0,3,0,0))

# plot conditional importance
resultsOutFluct.conimp <- arrange(resultsOutFluct$conditional_permimp, desc(row_number()))

pBar <- ggplot(resultsOutFluct.conimp) + geom_col(aes(x=factor(rownames(resultsOutFluct.conimp), levels=rownames(resultsOutFluct.conimp)), y=CondPermImp), fill=mycolours[4], width=0.5) + labs(x="Variable", y="Conditional variable permutation importance (mean change in annoyance)") + ggtitle("Fluctuation strength") + theme(text = element_text(family = "serif"), panel.grid=element_line(color = rgb(235, 235, 235, 100, maxColorValue = 255), linewidth = 0.25, linetype = 2)) + coord_flip()
pBar

if (saveplots){
  ggsave(filename="PtsABdAnnoyMnFluctConPermimp.svg", width=8, height=3.5, path=file.path(outFigPath, "svg"))
  unlink("PtsABdAnnoyMnFluctConPermimp.svg")
  
  ggsave(filename="PtsABdAnnoyMnFluctConPermimp.pdf", width=8, height=3.5, path=file.path(outFigPath, "pdf"))
  unlink("PtsABdAnnoyMnFluctConPermimp.pdf")
}

Selected metric

fluctVar <- "UASFluctECMA10ExBin"

Roughness

Set variables

# Roughness
iVars <- c(absVar, eventVar, ambVar, "UASRoughECMA10ExBin", "UASRoughECMA05ExBin", "UASRoughFZ10ExMaxLR", "UASRoughFZ05ExMaxLR", "UASRoughDW10ExMaxLR", "UASRoughDW05ExMaxLR")
dVar <- "dAnnoyMean"

seeds <- c(4701, 52187, 16589, 65217, 16893)

Hyperparameter tuning

p <- mtryTune(dataIn=stimDataNum, iVars=iVars, dVar=dVar, seeds=seeds[1:2],
             ntrees=ntrees, minsplit=minsplit, minbucket=minbucket)

p

Selected hyperparameters

ntree <- 1001
mtry <- as.integer(length(iVars)/1.8)

Run model

Train preliminary model

nperm <- 5

resultsOutRough <- crfReg(dataIn=stimDataNum, iVars=iVars, dVar=dVar, seeds=seeds[1:2], ntree=ntree, mtry=mtry, permImpCondThres=permImpCondThres, nperm=nperm, minsplit=minsplit, minbucket=minbucket)

# print model prediction results
resultsOutRough$OOB_RMSE

[1] 0.6741884

resultsOutRough$OOB_MAE

[1] 0.5361558

resultsOutRough$Rsquared

[1] 0.8539761

Train multiple seeds model

resultsOutRough <- multi_crfReg(dataIn=stimDataNum, iVars=iVars, dVar=dVar, seeds=seeds, ntree=ntree, mtry=mtry, permImpCondThres=permImpCondThres, nperm=nperm, minsplit=minsplit, minbucket=minbucket)

# print model prediction results
resultsOutRough$OOB_RMSE

[1] 0.6614793

resultsOutRough$OOB_MAE

[1] 0.5271245

resultsOutRough$Rsquared

[1] 0.8617422

# store results
resdAnnoyMnFitAB['Abs rough', 'RMSE'] <- resultsOutRough$OOB_RMSE
resdAnnoyMnFitAB['Abs rough', 'MAE'] <- resultsOutRough$OOB_MAE
resdAnnoyMnFitAB['Abs rough', 'Rsquared'] <- resultsOutRough$Rsquared
resdAnnoyMnPermImpAB$AbsRough <- resultsOutRough$conditional_permimp

Plot results

par(mai=c(0,3,0,0))

# plot conditional importance
resultsOutRough.conimp <- arrange(resultsOutRough$conditional_permimp, desc(row_number()))

pBar <- ggplot(resultsOutRough.conimp) + geom_col(aes(x=factor(rownames(resultsOutRough.conimp), levels=rownames(resultsOutRough.conimp)), y=CondPermImp), fill=mycolours[5], width=0.5) + labs(x="Variable", y="Conditional variable permutation importance (mean change in annoyance)") + ggtitle("Roughness") + theme(text = element_text(family = "serif"), panel.grid=element_line(color = rgb(235, 235, 235, 100, maxColorValue = 255), linewidth = 0.25, linetype = 2)) + coord_flip()
pBar

if (saveplots){
  ggsave(filename="PtsABdAnnoyMnRoughConPermimp.svg", width=8, height=2.9, path=file.path(outFigPath, "svg"))
  unlink("PtsABdAnnoyMnRoughConPermimp.svg")
  
  ggsave(filename="PtsABdAnnoyMnRoughConPermimp.pdf", width=8, height=2.9, path=file.path(outFigPath, "pdf"))
  unlink("PtsABdAnnoyMnRoughConPermimp.pdf")
}

Selected metric

roughVar <- "UASRoughFZ05ExMaxLR"

Impulsiveness

Set variables

# Impulsiveness
iVars <- c(absVar, eventVar, ambVar, "UASImpulsSHMAvgMaxLR", "UASImpulsSHM05ExMaxLR", "UASImpulsSHMPowAvgMaxLR", "UASImpulsLoudWZAvgMaxLR", "UASImpulsLoudWZ05ExMaxLR", "UASImpulsLoudWZPowAvgMaxLR", "UASImpulsLoudWECMAAvgBin", "UASImpulsLoudWECMA05ExBin", "UASImpulsLoudWECMAPowAvgBin")
dVar <- "dAnnoyMean"

seeds <- c(8495, 59867, 5416, 9843, 86)

Hyperparameter tuning

p <- mtryTune(dataIn=stimDataNum, iVars=iVars, dVar=dVar, seeds=seeds[1:2],
             ntrees=ntrees, minsplit=minsplit, minbucket=minbucket)

p

Selected hyperparameters

ntree <- 4001
mtry <- as.integer(length(iVars)/1.25)

Run model

Train preliminary model

nperm <- 5

resultsOutImpuls <- crfReg(dataIn=stimDataNum, iVars=iVars, dVar=dVar, seeds=seeds[1:2], ntree=ntree, mtry=mtry, permImpCondThres=permImpCondThres, nperm=nperm, minsplit=minsplit, minbucket=minbucket)

# print model prediction results
resultsOutImpuls$OOB_RMSE

[1] 0.6721437

resultsOutImpuls$OOB_MAE

[1] 0.5263928

resultsOutImpuls$Rsquared

[1] 0.8394367

Train multiple seeds model

resultsOutImpuls <- multi_crfReg(dataIn=stimDataNum, iVars=iVars, dVar=dVar, seeds=seeds, ntree=ntree, mtry=mtry, permImpCondThres=permImpCondThres, nperm=nperm, minsplit=minsplit, minbucket=minbucket)

# print model prediction results
resultsOutImpuls$OOB_RMSE

[1] 0.671666

resultsOutImpuls$OOB_MAE

[1] 0.525632

resultsOutImpuls$Rsquared

[1] 0.8396473

# store results
resdAnnoyMnFitAB['Abs impuls', 'RMSE'] <- resultsOutImpuls$OOB_RMSE
resdAnnoyMnFitAB['Abs impuls', 'MAE'] <- resultsOutImpuls$OOB_MAE
resdAnnoyMnFitAB['Abs impuls', 'Rsquared'] <- resultsOutImpuls$Rsquared
resdAnnoyMnPermImpAB$AbsImpuls <- resultsOutImpuls$conditional_permimp

Plot results

par(mai=c(0,3,0,0))

# plot conditional importance
resultsOutImpuls.conimp <- arrange(resultsOutImpuls$conditional_permimp, desc(row_number()))

pBar <- ggplot(resultsOutImpuls.conimp) + geom_col(aes(x=factor(rownames(resultsOutImpuls.conimp), levels=rownames(resultsOutImpuls.conimp)), y=CondPermImp), fill=mycolours[6], width=0.5) + labs(x="Variable", y="Conditional variable permutation importance (mean change in annoyance)") + ggtitle("Impulsiveness") + theme(text = element_text(family = "serif"), panel.grid=element_line(color = rgb(235, 235, 235, 100, maxColorValue = 255), linewidth = 0.25, linetype = 2)) + coord_flip()
pBar

if (saveplots){
  ggsave(filename="PtsABdAnnoyMnImpulsConPermimp.svg", width=8, height=3.8, path=file.path(outFigPath, "svg"))
  unlink("PtsABdAnnoyMnImpulsConPermimp.svg")
  
  ggsave(filename="PtsABdAnnoyMnImpulsConPermimp.pdf", width=8, height=3.8, path=file.path(outFigPath, "pdf"))
  unlink("PtsABdAnnoyMnImpulsConPermimp.pdf")
}

Selected metric

impulsVar <- "UASImpulsLoudWZAvgMaxLR"

SQM and loudness comparison

Now the highest importance SQMs are ranked against each other, controlling for UAS loudness and ambient LAeq.

Include tonal loudness

Set variables

iVars <- c(absVar, eventVar, ambVar, sharpVar, tonLdVar, fluctVar, roughVar, impulsVar)
dVar <- "dAnnoyMean"

seeds <- c(70498, 4, 14986, 453, 864)

Hyperparameter tuning

p <- mtryTune(dataIn=stimDataNum, iVars=iVars, dVar=dVar, seeds=seeds[1:2],
             ntrees=ntrees, minsplit=minsplit, minbucket=minbucket)

p

Selected hyperparameters

ntree <- 2501
mtry <- as.integer(length(iVars)/1.6)

Run model

Train preliminary model

nperm <- 5

resultsOutSQMs1 <- crfReg(dataIn=stimDataNum, iVars=iVars, dVar=dVar, seeds=seeds[1:2], ntree=ntree, mtry=mtry, permImpCondThres=permImpCondThres, nperm=nperm, minsplit=minsplit, minbucket=minbucket)

# print model prediction results
resultsOutSQMs1$OOB_RMSE

[1] 0.6690392

resultsOutSQMs1$OOB_MAE

[1] 0.5373742

resultsOutSQMs1$Rsquared

[1] 0.8465967

Train multiple seeds model

resultsOutSQMs1 <- multi_crfReg(dataIn=stimDataNum, iVars=iVars, dVar=dVar, seeds=seeds, ntree=ntree, mtry=mtry, permImpCondThres=permImpCondThres, nperm=nperm, minsplit=minsplit, minbucket=minbucket)

# print model prediction results
resultsOutSQMs1$OOB_RMSE

[1] 0.6651781

resultsOutSQMs1$OOB_MAE

[1] 0.5349845

resultsOutSQMs1$Rsquared

[1] 0.8487605

# store results
resdAnnoyMnFitAB['Abs SQMs inc tonal loud', 'RMSE'] <- resultsOutSQMs1$OOB_RMSE
resdAnnoyMnFitAB['Abs SQMs inc tonal loud', 'MAE'] <- resultsOutSQMs1$OOB_MAE
resdAnnoyMnFitAB['Abs SQMs inc tonal loud', 'Rsquared'] <- resultsOutSQMs1$Rsquared
resdAnnoyMnPermImpAB$AbsSQMs1 <- resultsOutSQMs1$conditional_permimp

Plot results

par(mai=c(0,3,0,0))

# plot conditional importance
resultsOutSQMs1.conimp <- arrange(resultsOutSQMs1$conditional_permimp, desc(row_number()))

pBar <- ggplot(resultsOutSQMs1.conimp) + geom_col(aes(x=factor(rownames(resultsOutSQMs1.conimp), levels=rownames(resultsOutSQMs1.conimp)), y=CondPermImp), fill=mycolours[7], width=0.5) + labs(x="Variable", y="Conditional variable permutation importance (mean change in annoyance)") + theme(text = element_text(family = "serif"), panel.grid=element_line(color = rgb(235, 235, 235, 100, maxColorValue = 255), linewidth = 0.25, linetype = 2)) + coord_flip(ylim=c(0, 1.3))
pBar

if (saveplots){
  ggsave(filename="PtsABdAnnoyMnAbsSQMsTonLdConPermimp.svg", width=8, height=2.4, path=file.path(outFigPath, "svg"))
  unlink("PtsABdAnnoyMnAbsSQMsTonLdConPermimp.svg")
  
  ggsave(filename="PtsABdAnnoyMnAbsSQMsTonLdConPermimp.pdf", width=8, height=2.4, path=file.path(outFigPath, "pdf"))
  unlink("PtsABdAnnoyMnAbsSQMsTonLdConPermimp.pdf")
}

Exclude tonal loudness

Set variables

iVars <- c(absVar, eventVar, ambVar, sharpVar, tonalVar, fluctVar, roughVar, impulsVar)
dVar <- "dAnnoyMean"

seeds <- c(546, 57203, 270835, 60592, 8094)

Hyperparameter tuning

p <- mtryTune(dataIn=stimDataNum, iVars=iVars, dVar=dVar, seeds=seeds[1:2],
             ntrees=ntrees, minsplit=minsplit, minbucket=minbucket)

p

Selected hyperparameters

ntree <- 1501
mtry <- as.integer(length(iVars)/1.6)

Run model

Train preliminary model

nperm <- 5

resultsOutSQMs2 <- crfReg(dataIn=stimDataNum, iVars=iVars, dVar=dVar, seeds=seeds[1:2], ntree=ntree, mtry=mtry, permImpCondThres=permImpCondThres, nperm=nperm, minsplit=minsplit, minbucket=minbucket)

# print model prediction results
resultsOutSQMs2$OOB_RMSE

[1] 0.6741422

resultsOutSQMs2$OOB_MAE

[1] 0.5413668

resultsOutSQMs2$Rsquared

[1] 0.8439303

Train multiple seeds model

resultsOutSQMs2 <- multi_crfReg(dataIn=stimDataNum, iVars=iVars, dVar=dVar, seeds=seeds, ntree=ntree, mtry=mtry, permImpCondThres=permImpCondThres, nperm=nperm, minsplit=minsplit, minbucket=minbucket)

# print model prediction results
resultsOutSQMs2$OOB_RMSE

[1] 0.6606138

resultsOutSQMs2$OOB_MAE

[1] 0.5318176

resultsOutSQMs2$Rsquared

[1] 0.8522687

# store results
resdAnnoyMnFitAB['Abs SQMs no tonal loud', 'RMSE'] <- resultsOutSQMs2$OOB_RMSE
resdAnnoyMnFitAB['Abs SQMs no tonal loud', 'MAE'] <- resultsOutSQMs2$OOB_MAE
resdAnnoyMnFitAB['Abs SQMs no tonal loud', 'Rsquared'] <- resultsOutSQMs2$Rsquared
resdAnnoyMnPermImpAB$AbsSQMs2 <- resultsOutSQMs2$conditional_permimp

Plot results

par(mai=c(0,3,0,0))

# plot conditional importance
resultsOutSQMs2.conimp <- arrange(resultsOutSQMs2$conditional_permimp, desc(row_number()))

pBar <- ggplot(resultsOutSQMs2.conimp) + geom_col(aes(x=factor(rownames(resultsOutSQMs2.conimp), levels=rownames(resultsOutSQMs2.conimp)), y=CondPermImp), fill=mycolours[7], width=0.5) + labs(x="Variable", y="Conditional variable permutation importance (mean change in annoyance)") + theme(text = element_text(family = "serif"), panel.grid=element_line(color = rgb(235, 235, 235, 100, maxColorValue = 255), linewidth = 0.25, linetype = 2)) + coord_flip(ylim=c(0, 1.3))
pBar

if (saveplots){
  ggsave(filename="PtsABdAnnoyMnAbsSQMsNoTonLdConPermimp.svg", width=8, height=2.4, path=file.path(outFigPath, "svg"))
  unlink("PtsABdAnnoyMnAbsSQMsNoTonLdConPermimp.svg")
  
  ggsave(filename="PtsABdAnnoyMnAbsSQMsNoTonLdConPermimp.pdf", width=8, height=2.4, path=file.path(outFigPath, "pdf"))
  unlink("PtsABdAnnoyMnAbsSQMsNoTonLdConPermimp.pdf")
}

Psychoacoustic annoyance metrics

Set variables

iVars <- c(ambVar, "UASPsychAnnoyWidmann", "UASPsychAnnoyMore", "UASPsychAnnoyDi", "UASPsychAnnoyTorija", "UASPsychAnnoyWillemsen", "UASPsychAnnoyBoucher")
dVar <- "dAnnoyMean"

seeds <- c(829, 9, 190, 4564, 924707824)

Hyperparameter tuning

p <- mtryTune(dataIn=stimDataNum, iVars=iVars, dVar=dVar, seeds=seeds[1:2],
             ntrees=ntrees, minsplit=minsplit, minbucket=minbucket)

p

Selected hyperparameters

ntree <- 1501
mtry <- as.integer(length(iVars)/1.75)

Run model

Train preliminary model

nperm <- 5

resultsOutPA <- crfReg(dataIn=stimDataNum, iVars=iVars, dVar=dVar, seeds=seeds[1:2], ntree=ntree, mtry=mtry, permImpCondThres=permImpCondThres, nperm=nperm, minsplit=minsplit, minbucket=minbucket)

# print model prediction results
resultsOutPA$OOB_RMSE

[1] 0.6959625

resultsOutPA$OOB_MAE

[1] 0.5694917

resultsOutPA$Rsquared

[1] 0.8202226

Train multiple seeds model

resultsOutPA <- multi_crfReg(dataIn=stimDataNum, iVars=iVars, dVar=dVar, seeds=seeds, ntree=ntree, mtry=mtry, permImpCondThres=permImpCondThres, nperm=nperm, minsplit=minsplit, minbucket=minbucket)

# print model prediction results
resultsOutPA$OOB_RMSE

[1] 0.6987601

resultsOutPA$OOB_MAE

[1] 0.5691674

resultsOutPA$Rsquared

[1] 0.8188134

# store results
resdAnnoyMnFitAB['Psychoacoustic annoyance', 'RMSE'] <- resultsOutPA$OOB_RMSE
resdAnnoyMnFitAB['Psychoacoustic annoyance', 'MAE'] <- resultsOutPA$OOB_MAE
resdAnnoyMnFitAB['Psychoacoustic annoyance', 'Rsquared'] <- resultsOutPA$Rsquared
resdAnnoyMnPermImpAB$AbsPA <- resultsOutPA$conditional_permimp

Plot results

par(mai=c(0,3,0,0))

# plot conditional importance
resultsOutPA.conimp <- arrange(resultsOutPA$conditional_permimp, desc(row_number()))

pBar <- ggplot(resultsOutPA.conimp) + geom_col(aes(x=factor(rownames(resultsOutPA.conimp), levels=rownames(resultsOutPA.conimp)), y=CondPermImp), fill=mycolours[10], width=0.5) + labs(x="Variable", y="Conditional variable permutation importance (mean change in annoyance)") + theme(text = element_text(family = "serif"), panel.grid=element_line(color = rgb(235, 235, 235, 100, maxColorValue = 255), linewidth = 0.25, linetype = 2)) + coord_flip(ylim=c(0, 1.8))
pBar

if (saveplots){
  ggsave(filename="PtsABdAnnoyMnAbsPAConPermimp.svg", width=8, height=2.4, path=file.path(outFigPath, "svg"))
  unlink("PtsABdAnnoyMnAbsPAConPermimp.svg")
  
  ggsave(filename="PtsABdAnnoyMnAbsPAConPermimp.pdf", width=8, height=2.4, path=file.path(outFigPath, "pdf"))
  unlink("PtsABdAnnoyMnAbsPAConPermimp.pdf")
}

Difference variables

Next, the difference metrics are analysed,

Set variables

iVars <- c(names(stimDataNum)[which(colnames(stimDataNum)=="LAeqLAF90diff"):
                               which(colnames(stimDataNum)=="dPsychAnnoyBoucher")],
           "UASEvents", "SNRlevel", "IntermitRatioC2MaxLR", "IntermitRatioC3MaxLR", "IntermitRatioC5MaxLR")
dVar <- "dAnnoyMean"

seeds <- c(568392, 498, 4089, 78132, 741809)

Hyperparameter tuning

p <- mtryTune(dataIn=stimDataNum, iVars=iVars, dVar=dVar, seeds=seeds[1:2],
             ntrees=ntrees, minsplit=minsplit, minbucket=minbucket)

p

Selected hyperparameters

ntree <- 4001
mtry <- as.integer(length(iVars)/2.25)

Run model

Train preliminary model

nperm <- 5

resultsOutDiffs <- crfReg(dataIn=stimDataNum, iVars=iVars, dVar=dVar, seeds=seeds[1:2], ntree=ntree, mtry=mtry, permImpCondThres=permImpCondThres, nperm=nperm, minsplit=minsplit, minbucket=minbucket)

# print model prediction results
resultsOutDiffs$OOB_RMSE

[1] 0.5275899

resultsOutDiffs$OOB_MAE

[1] 0.4075381

resultsOutDiffs$Rsquared

[1] 0.8912516

Train multiple seeds model

resultsOutDiffs <- multi_crfReg(dataIn=stimDataNum, iVars=iVars, dVar=dVar, seeds=seeds, ntree=ntree, mtry=mtry, permImpCondThres=permImpCondThres, nperm=nperm, minsplit=minsplit, minbucket=minbucket)

# print model prediction results
resultsOutDiffs$OOB_RMSE

[1] 0.5265303

resultsOutDiffs$OOB_MAE

[1] 0.4070696

resultsOutDiffs$Rsquared

[1] 0.8916423

# store results
resdAnnoyMnFitAB['Diff vars', 'RMSE'] <- resultsOutDiffs$OOB_RMSE
resdAnnoyMnFitAB['Diff vars', 'MAE'] <- resultsOutDiffs$OOB_MAE
resdAnnoyMnFitAB['Diff vars', 'Rsquared'] <- resultsOutDiffs$Rsquared
resdAnnoyMnPermImpAB$DiffVars <- resultsOutDiffs$conditional_permimp

Plot results

par(mai=c(0,3,0,0))

# plot conditional importance
resultsOutDiffs.conimp <- arrange(resultsOutDiffs$conditional_permimp, desc(row_number()))

pBar <- ggplot(resultsOutDiffs.conimp) + geom_col(aes(x=factor(rownames(resultsOutDiffs.conimp), levels=rownames(resultsOutDiffs.conimp)), y=CondPermImp), fill=mycolours[8], width=0.5) + labs(x="Variable", y="Conditional variable permutation importance (mean change in annoyance)") + theme(text = element_text(family = "serif"), panel.grid=element_line(color = rgb(235, 235, 235, 100, maxColorValue = 255), linewidth = 0.25, linetype = 2)) + coord_flip()
pBar

if (saveplots){
  ggsave(filename="PtsABdAnnoyMnDiffVarsConPermimp.svg", width=8, height=18, path=file.path(outFigPath, "svg"))
  unlink("PtsABdAnnoyMnDiffVarsConPermimp.svg")

  ggsave(filename="PtsABdAnnoyMnDiffVarsConPermimp.pdf", width=8, height=18, path=file.path(outFigPath, "pdf"))
  unlink("PtsABdAnnoyMnDiffVarsConPermimp.pdf")
}

# Plot only positive values
resultsOutDiffs.conimpPtv <- resultsOutDiffs.conimp |>
                                        rownames_to_column('Metric') |>
                                            filter_if(is.numeric, all_vars(. > 0)) |>
                                                  column_to_rownames('Metric')

pBar <- ggplot(resultsOutDiffs.conimpPtv) + geom_col(aes(x=factor(rownames(resultsOutDiffs.conimpPtv), levels=rownames(resultsOutDiffs.conimpPtv)), y=CondPermImp), fill=mycolours[8], width=0.5) + labs(x="Variable", y="Conditional variable permutation importance (mean change in annoyance)") + theme(text = element_text(family = "serif"), panel.grid=element_line(color = rgb(235, 235, 235, 100, maxColorValue = 255), linewidth = 0.25, linetype = 2)) + coord_flip()
pBar

if (saveplots){
  ggsave(filename="PtsABdAnnoyMnDiffVarsConPermimpPtv.svg", width=8, height=12, path=file.path(outFigPath, "svg"))
  unlink("PtsABdAnnoyMnDiffVarsConPermimpPtv.svg")

  ggsave(filename="PtsABdAnnoyMnDiffVarsConPermimpPtv.pdf", width=8, height=12, path=file.path(outFigPath, "pdf"))
  unlink("PtsABdAnnoyMnDiffVarsConPermimpPtv.pdf")
}

# Plot only values within 1% of the maximum
resultsOutDiffs.conimp1pc <- resultsOutDiffs.conimp |>
                                        rownames_to_column('Metric') |>
                                            filter_if(is.numeric, all_vars(. > max(resultsOutDiffs.conimp)/100)) |>
                                                  column_to_rownames('Metric')

pBar <- ggplot(resultsOutDiffs.conimp1pc) + geom_col(aes(x=factor(rownames(resultsOutDiffs.conimp1pc), levels=rownames(resultsOutDiffs.conimp1pc)), y=CondPermImp), fill=mycolours[8], width=0.5) + labs(x="Variable", y="Conditional variable permutation importance (mean change in annoyance)") + theme(text = element_text(family = "serif"), panel.grid=element_line(color = rgb(235, 235, 235, 100, maxColorValue = 255), linewidth = 0.25, linetype = 2)) + coord_flip()
pBar

if (saveplots){
  ggsave(filename="PtsABdAnnoyMnDiffVarsConPermimp1pc.svg", width=8, height=5, path=file.path(outFigPath, "svg"))
  unlink("PtsABdAnnoyMnDiffVarsConPermimp1pc.svg")

  ggsave(filename="PtsABdAnnoyMnDiffVarsConPermimp1pc.pdf", width=8, height=5, path=file.path(outFigPath, "pdf"))
  unlink("PtsABdAnnoyMnDiffVarsConPermimp1pc.pdf")
}

# Plot only values within 1% of the maximum
resultsOutDiffs.conimp1pc <- resultsOutDiffs.conimp |>
                                        rownames_to_column('Metric') |>
                                            filter_if(is.numeric, all_vars(. > max(resultsOutDiffs.conimp)/100)) |>
                                                  column_to_rownames('Metric')

pBar <- ggplot(resultsOutDiffs.conimp1pc) + geom_col(aes(x=factor(rownames(resultsOutDiffs.conimp1pc), levels=rownames(resultsOutDiffs.conimp1pc)), y=CondPermImp), fill=mycolours[8], width=0.5) + labs(x="Variable", y="Conditional variable permutation\nimportance (mean change in annoyance)") + theme(text = element_text(family = "serif"), axis.title.x=element_text(family = "serif", size=9), panel.grid=element_line(color = rgb(235, 235, 235, 100, maxColorValue = 255), linewidth = 0.25, linetype = 2)) + coord_flip()
pBar

filename = "PtsABdAnnoyMnDiffVarsConPermimp1pcNw"

if (saveplots){
  ggsave(filename=paste(filename, ".svg"), width=4, height=5, path=file.path(outFigPath, "svg"))
  unlink(paste(filename, ".svg"))

  ggsave(filename=paste(filename, ".pdf"), width=4, height=5, path=file.path(outFigPath, "pdf"))
  unlink(paste(filename, ".svg"))
}

Selected metric

diffVar <- "LAELAF50diff"

dSQM analysis

Individual SQMs

dSharpness

Set variables

iVars <- c(diffVar, eventVar, "dSharpAurISO3PowAvgBin", "dSharpAurISO305ExBin", "dSharpAurSHMPowAvgBin", "dSharpAurSHM05ExBin", "dTonShpAurSHMPowAvgBin", "dTonShpAurSHM05ExBin", "PartTonShpAurSHMPowAvgBin",
           "PartTonShpAurSHM05ExBin")
dVar <- "dAnnoyMean"

seeds <- c(84194, 905, 64815, 928054, 625091)

Hyperparameter tuning

p <- mtryTune(dataIn=stimDataNum, iVars=iVars, dVar=dVar, seeds=seeds[1:2],
             ntrees=ntrees, minsplit=minsplit, minbucket=minbucket)

p

Selected hyperparameters

ntree <-5501
mtry <- as.integer(length(iVars)/1.25)

Run model

Train preliminary model

nperm <- 5

resultsOutSharp <- crfReg(dataIn=stimDataNum, iVars=iVars, dVar=dVar, seeds=seeds[1:2], ntree=ntree, mtry=mtry, permImpCondThres=permImpCondThres, nperm=nperm, minsplit=minsplit, minbucket=minbucket)

# print model prediction results
resultsOutSharp$OOB_RMSE

[1] 0.5364699

resultsOutSharp$OOB_MAE

[1] 0.4069908

resultsOutSharp$Rsquared

[1] 0.8883255

Train multiple seeds model

resultsOutSharp <- multi_crfReg(dataIn=stimDataNum, iVars=iVars, dVar=dVar, seeds=seeds, ntree=ntree, mtry=mtry, permImpCondThres=permImpCondThres, nperm=nperm, minsplit=minsplit, minbucket=minbucket)

# print model prediction results
resultsOutSharp$OOB_RMSE

[1] 0.53737

resultsOutSharp$OOB_MAE

[1] 0.4076505

resultsOutSharp$Rsquared

[1] 0.8880165

# store results
resdAnnoyMnFitAB['Diff sharp', 'RMSE'] <- resultsOutSharp$OOB_RMSE
resdAnnoyMnFitAB['Diff sharp', 'MAE'] <- resultsOutSharp$OOB_MAE
resdAnnoyMnFitAB['Diff sharp', 'Rsquared'] <- resultsOutSharp$Rsquared
resdAnnoyMnPermImpAB$DiffSharp <- resultsOutSharp$conditional_permimp

Plot results

par(mai=c(0,3,0,0))

# plot conditional importance
resultsOutSharp.conimp <- arrange(resultsOutSharp$conditional_permimp, desc(row_number()))

pBar <- ggplot(resultsOutSharp.conimp) + geom_col(aes(x=factor(rownames(resultsOutSharp.conimp), levels=rownames(resultsOutSharp.conimp)), y=CondPermImp), fill=mycolours[2], width=0.5) + labs(x="Variable", y="Conditional variable permutation importance (mean change in annoyance)") + ggtitle("dSharpness") + theme(text = element_text(family = "serif"), panel.grid=element_line(color = rgb(235, 235, 235, 100, maxColorValue = 255), linewidth = 0.25, linetype = 2)) + coord_flip()
pBar

if (saveplots){
  ggsave(filename="PtsABdAnnoyMndSharpConPermimp.svg", width=8, height=3.4, path=file.path(outFigPath, "svg"))
  unlink("PtsABdAnnoyMndSharpConPermimp.svg")
  
  ggsave(filename="PtsABdAnnoyMndSharpConPermimp.pdf", width=8, height=3.4, path=file.path(outFigPath, "pdf"))
  unlink("PtsABdAnnoyMndSharpConPermimp.pdf")
}

Selected metric

dSharpVar <- "dSharpAurISO3PowAvgBin"

dTonal loudness and dtonality

Set variables

iVars <- c(diffVar, eventVar, "dTonalECMAAvgMaxLR", "dTonalSHMInt05ExMaxLR", "dTonalSHMIntAvgMaxLR", "dTonalECMA05ExMaxLR", "dTonalAwSHMAvgMaxLR",  "dTonalAwSHM05ExMaxLR", "dTonalAwSHMIntAvgMaxLR",   "dTonalAwSHMInt05ExMaxLR", "dTonLdECMAPowAvgBin", "dTonLdECMA05ExBin", "dTonShpAurSHMPowAvgBin",
           "dTonShpAurSHM05ExBin", "PartTonLdSHMPowAvgBin")
dVar <- "dAnnoyMean"

seeds <- c(561684, 104798, 1536, 48, 48561)

Hyperparameter tuning

p <- mtryTune(dataIn=stimDataNum, iVars=iVars, dVar=dVar, seeds=seeds[1:2],
             ntrees=ntrees, minsplit=minsplit, minbucket=minbucket)

p

Selected hyperparameters

ntree <- 1001
mtry <- as.integer(length(iVars)/1.75)

Run model

Train preliminary model

# Tonality with tonal loudness

nperm <- 5

resultsOutTonal1 <- crfReg(dataIn=stimDataNum, iVars=iVars, dVar=dVar, seeds=seeds[1:2], ntree=ntree, mtry=mtry, permImpCondThres=permImpCondThres, nperm=nperm, minsplit=minsplit, minbucket=minbucket)

# print model prediction results
resultsOutTonal1$OOB_RMSE

[1] 0.531723

resultsOutTonal1$OOB_MAE

[1] 0.401651

resultsOutTonal1$Rsquared

[1] 0.8899811

Train multiple seeds model

# Tonality with tonal loudness

resultsOutTonal1 <- multi_crfReg(dataIn=stimDataNum, iVars=iVars, dVar=dVar, seeds=seeds, ntree=ntree, mtry=mtry, permImpCondThres=permImpCondThres, nperm=nperm, minsplit=minsplit, minbucket=minbucket)

# print model prediction results
resultsOutTonal1$OOB_RMSE

[1] 0.533725

resultsOutTonal1$OOB_MAE

[1] 0.4030735

resultsOutTonal1$Rsquared

[1] 0.889263

# store results
resdAnnoyMnFitAB['Diff tonal inc loud', 'RMSE'] <- resultsOutTonal1$OOB_RMSE
resdAnnoyMnFitAB['Diff tonal inc loud', 'MAE'] <- resultsOutTonal1$OOB_MAE
resdAnnoyMnFitAB['Diff tonal inc loud', 'Rsquared'] <- resultsOutTonal1$Rsquared
resdAnnoyMnPermImpAB$DiffTonal1 <- resultsOutTonal1$conditional_permimp

Plot results

par(mai=c(0,3,0,0))

# plot conditional importance
resultsOutTonal1.conimp <- arrange(resultsOutTonal1$conditional_permimp, desc(row_number()))

pBar <- ggplot(resultsOutTonal1.conimp) + geom_col(aes(x=factor(rownames(resultsOutTonal1.conimp), levels=rownames(resultsOutTonal1.conimp)), y=CondPermImp), fill=mycolours[3], width=0.5) + labs(x="Variable", y="Conditional variable permutation importance (mean change in annoyance)") + ggtitle("dTonality inc. dtonal loudness") + theme(text = element_text(family = "serif"), panel.grid=element_line(color = rgb(235, 235, 235, 100, maxColorValue = 255), linewidth = 0.25, linetype = 2)) + coord_flip(ylim=c(0, 2.2))
pBar

if (saveplots){
  ggsave(filename="PtsABdAnnoyMndTonalLdConPermimp.svg", width=8, height=3.6, path=file.path(outFigPath, "svg"))
  unlink("PtsABdAnnoyMndTonalLdConPermimp.svg")
  
  ggsave(filename="PtsABdAnnoyMndTonalLdConPermimp.pdf", width=8, height=3.6, path=file.path(outFigPath, "pdf"))
  unlink("PtsABdAnnoyMndTonalLdConPermimp.pdf")
}

Selected metric

dTonLdVar <- "dTonLdECMAPowAvgBin"

dTonality without dtonal loudness

Set variables

iVars <- c(diffVar, eventVar, "dTonalECMAAvgMaxLR", "dTonalSHMInt05ExMaxLR", "dTonalSHMIntAvgMaxLR", "dTonalECMA05ExMaxLR", "dTonalAwSHMAvgMaxLR",  "dTonalAwSHM05ExMaxLR", "dTonalAwSHMIntAvgMaxLR",   "dTonalAwSHMInt05ExMaxLR")
dVar <- "dAnnoyMean"

seeds <- c(410865, 2954, 70812, 203, 7984)

Hyperparameter tuning

p <- mtryTune(dataIn=stimDataNum, iVars=iVars, dVar=dVar, seeds=seeds[1:2],
             ntrees=ntrees, minsplit=minsplit, minbucket=minbucket)

p

Selected hyperparameters

ntree <- 5501
mtry <- as.integer(length(iVars)/1.25)

Run model

Train preliminary model

# Tonality

nperm <- 5

resultsOutTonal2 <- crfReg(dataIn=stimDataNum, iVars=iVars, dVar=dVar, seeds=seeds[1:2],
                           ntree=ntree, mtry=mtry, permImpCondThres=permImpCondThres, nperm=nperm, minsplit=minsplit, minbucket=minbucket)

# print model prediction results
resultsOutTonal2$OOB_RMSE

[1] 0.5363242

resultsOutTonal2$OOB_MAE

[1] 0.409025

resultsOutTonal2$Rsquared

[1] 0.8878642

Train multiple seeds model

# Tonality

resultsOutTonal2 <- multi_crfReg(dataIn=stimDataNum, iVars=iVars, dVar=dVar, seeds=seeds, ntree=ntree, mtry=mtry, permImpCondThres=permImpCondThres, nperm=nperm, minsplit=minsplit, minbucket=minbucket)

# print model prediction results
resultsOutTonal2$OOB_RMSE

[1] 0.5352787

resultsOutTonal2$OOB_MAE

[1] 0.4084578

resultsOutTonal2$Rsquared

[1] 0.8883184

# store results
resdAnnoyMnFitAB['Diff tonal no loud', 'RMSE'] <- resultsOutTonal2$OOB_RMSE
resdAnnoyMnFitAB['Diff tonal no loud', 'MAE'] <- resultsOutTonal2$OOB_MAE
resdAnnoyMnFitAB['Diff tonal no loud', 'Rsquared'] <- resultsOutTonal2$Rsquared
resdAnnoyMnPermImpAB$DiffTonal2 <- resultsOutTonal2$conditional_permimp

Plot results

par(mai=c(0,3,0,0))

# plot conditional importance
resultsOutTonal2.conimp <- arrange(resultsOutTonal2$conditional_permimp, desc(row_number()))

pBar <- ggplot(resultsOutTonal2.conimp) + geom_col(aes(x=factor(rownames(resultsOutTonal2.conimp), levels=rownames(resultsOutTonal2.conimp)), y=CondPermImp), fill=mycolours[3], width=0.5) + labs(x="Variable", y="Conditional variable permutation importance (mean change in annoyance)") + ggtitle("dTonality w/o tonal loudness") + theme(text = element_text(family = "serif"), panel.grid=element_line(color = rgb(235, 235, 235, 100, maxColorValue = 255), linewidth = 0.25, linetype = 2)) + coord_flip(ylim=c(0, 2.2))
pBar

if (saveplots){
  ggsave(filename="PtsABdAnnoyMndTonalConPermimp.svg", width=8, height=2.8, path=file.path(outFigPath, "svg"))
  unlink("PtsABdAnnoyMndTonalConPermimp.svg")
  
  ggsave(filename="PtsABdAnnoyMndTonalConPermimp.pdf", width=8, height=2.8, path=file.path(outFigPath, "pdf"))
  unlink("PtsABdAnnoyMndTonalConPermimp.pdf")
}

Selected metric

dTonalVar <- "dTonalSHMIntAvgMaxLR"

dFluctuation strength

Set variables

# Fluctuation strength
iVars <- c(diffVar, eventVar, "dFluctECMA10ExBin", "dFluctECMA05ExBin", "dFluctOV10ExMaxLR", "dFluctOV05ExMaxLR")
dVar <- "dAnnoyMean"

seeds <- c(418657, 84, 1630, 18659, 3687)

Hyperparameter tuning

p <- mtryTune(dataIn=stimDataNum, iVars=iVars, dVar=dVar, seeds=seeds[1:2],
             ntrees=ntrees, minsplit=minsplit, minbucket=minbucket)

p

Selected hyperparameters

ntree <- 251
mtry <- as.integer(length(iVars)/2)

Run model

Train preliminary model

nperm <- 5

resultsOutFluct <- crfReg(dataIn=stimDataNum, iVars=iVars, dVar=dVar, seeds=seeds[1:2], ntree=ntree, mtry=mtry, permImpCondThres=permImpCondThres, nperm=nperm, minsplit=minsplit, minbucket=minbucket)

# print model prediction results
resultsOutFluct$OOB_RMSE

[1] 0.5306371

resultsOutFluct$OOB_MAE

[1] 0.4069944

resultsOutFluct$Rsquared

[1] 0.8917573

Train multiple seeds model

resultsOutFluct <- multi_crfReg(dataIn=stimDataNum, iVars=iVars, dVar=dVar, seeds=seeds, ntree=ntree, mtry=mtry, permImpCondThres=permImpCondThres, nperm=nperm, minsplit=minsplit, minbucket=minbucket)

# print model prediction results
resultsOutFluct$OOB_RMSE

[1] 0.5370137

resultsOutFluct$OOB_MAE

[1] 0.4133121

resultsOutFluct$Rsquared

[1] 0.8889464

# store results
resdAnnoyMnFitAB['Diff fluct', 'RMSE'] <- resultsOutFluct$OOB_RMSE
resdAnnoyMnFitAB['Diff fluct', 'MAE'] <- resultsOutFluct$OOB_MAE
resdAnnoyMnFitAB['Diff fluct', 'Rsquared'] <- resultsOutFluct$Rsquared
resdAnnoyMnPermImpAB$DiffFluct <- resultsOutFluct$conditional_permimp

Plot results

par(mai=c(0,3,0,0))

# plot conditional importance
resultsOutFluct.conimp <- arrange(resultsOutFluct$conditional_permimp, desc(row_number()))

pBar <- ggplot(resultsOutFluct.conimp) + geom_col(aes(x=factor(rownames(resultsOutFluct.conimp), levels=rownames(resultsOutFluct.conimp)), y=CondPermImp), fill=mycolours[4], width=0.5) + labs(x="Variable", y="Conditional variable permutation importance (mean change in annoyance)") + ggtitle("dFluctuation strength") + theme(text = element_text(family = "serif"), panel.grid=element_line(color = rgb(235, 235, 235, 100, maxColorValue = 255), linewidth = 0.25, linetype = 2)) + coord_flip()
pBar

if (saveplots){
  ggsave(filename="PtsABdAnnoyMndFluctConPermimp.svg", width=8, height=2, path=file.path(outFigPath, "svg"))
  unlink("PtsABdAnnoyMnFluctConPermimp.svg")
  
  ggsave(filename="PtsABdAnnoyMndFluctConPermimp.pdf", width=8, height=2, path=file.path(outFigPath, "pdf"))
  unlink("PtsABdAnnoyMnFluctConPermimp.pdf")
}

Selected metric

dFluctVar <- "dFluctECMA10ExBin"

dRoughness

Set variables

# Roughness
iVars <- c(diffVar, eventVar, "dRoughECMA10ExBin", "dRoughECMA05ExBin", "dRoughFZ10ExMaxLR", "dRoughFZ05ExMaxLR")
dVar <- "dAnnoyMean"

seeds <- c(69851, 85109, 410986, 1563, 896)

Hyperparameter tuning

p <- mtryTune(dataIn=stimDataNum, iVars=iVars, dVar=dVar, seeds=seeds[1:2],
             ntrees=ntrees, minsplit=minsplit, minbucket=minbucket)

p

Selected hyperparameters

ntree <- 5501
mtry <- as.integer(length(iVars)/2)

Run model

Train preliminary model

nperm <- 5

resultsOutRough <- crfReg(dataIn=stimDataNum, iVars=iVars, dVar=dVar, seeds=seeds[1:2], ntree=ntree, mtry=mtry, permImpCondThres=permImpCondThres, nperm=nperm, minsplit=minsplit, minbucket=minbucket)

# print model prediction results
resultsOutRough$OOB_RMSE

[1] 0.5909921

resultsOutRough$OOB_MAE

[1] 0.4655595

resultsOutRough$Rsquared

[1] 0.8792098

Train multiple seeds model

resultsOutRough <- multi_crfReg(dataIn=stimDataNum, iVars=iVars, dVar=dVar, seeds=seeds, ntree=ntree, mtry=mtry, permImpCondThres=permImpCondThres, nperm=nperm, minsplit=minsplit, minbucket=minbucket)

# print model prediction results
resultsOutRough$OOB_RMSE

[1] 0.5925856

resultsOutRough$OOB_MAE

[1] 0.4657394

resultsOutRough$Rsquared

[1] 0.8786725

# store results
resdAnnoyMnFitAB['Diff rough', 'RMSE'] <- resultsOutRough$OOB_RMSE
resdAnnoyMnFitAB['Diff rough', 'MAE'] <- resultsOutRough$OOB_MAE
resdAnnoyMnFitAB['Diff rough', 'Rsquared'] <- resultsOutRough$Rsquared
resdAnnoyMnPermImpAB$DiffRough <- resultsOutRough$conditional_permimp

Plot results

par(mai=c(0,3,0,0))

# plot conditional importance
resultsOutRough.conimp <- arrange(resultsOutRough$conditional_permimp, desc(row_number()))

pBar <- ggplot(resultsOutRough.conimp) + geom_col(aes(x=factor(rownames(resultsOutRough.conimp), levels=rownames(resultsOutRough.conimp)), y=CondPermImp), fill=mycolours[5], width=0.5) + labs(x="Variable", y="Conditional variable permutation importance (mean change in annoyance)") + ggtitle("dRoughness") + theme(text = element_text(family = "serif"), panel.grid=element_line(color = rgb(235, 235, 235, 100, maxColorValue = 255), linewidth = 0.25, linetype = 2)) + coord_flip()
pBar

if (saveplots){
  ggsave(filename="PtsABdAnnoyMndRoughConPermimp.svg", width=8, height=2, path=file.path(outFigPath, "svg"))
  unlink("PtsABdAnnoyMndRoughConPermimp.svg")
  
  ggsave(filename="PtsABdAnnoyMndRoughConPermimp.pdf", width=8, height=2, path=file.path(outFigPath, "pdf"))
  unlink("PtsABdAnnoyMndRoughConPermimp.pdf")
}

Selected metric

dRoughVar <- "dRoughFZ05ExMaxLR"

dImpulsiveness

Set variables

# Impulsiveness
iVars <- c(diffVar, eventVar, "dImpulsSHMAvgMaxLR", "dImpulsSHM05ExMaxLR", "dImpulsSHMPowAvgMaxLR",
           "dImpulsLoudWZAvgMaxLR", "dImpulsLoudWZ05ExMaxLR", "dImpulsLoudWZPowAvgMaxLR",
           "dImpulsLoudWECMAAvgBin", "dImpulsLoudWECMA05ExBin", "dImpulsLoudWECMAPowAvgBin")
dVar <- "dAnnoyMean"

seeds <- c(418659, 7805, 38475, 65834, 1653)

Hyperparameter tuning

p <- mtryTune(dataIn=stimDataNum, iVars=iVars, dVar=dVar, seeds=seeds[1:2],
             ntrees=ntrees, minsplit=minsplit, minbucket=minbucket)

p

Selected hyperparameters

ntree <- 1001
mtry <- as.integer(length(iVars)/1.25)

Run model

Train preliminary model

nperm <- 5

resultsOutImpuls <- crfReg(dataIn=stimDataNum, iVars=iVars, dVar=dVar, seeds=seeds[1:2], ntree=ntree, mtry=mtry, permImpCondThres=permImpCondThres, nperm=nperm, minsplit=minsplit, minbucket=minbucket)

# print model prediction results
resultsOutImpuls$OOB_RMSE

[1] 0.5340323

resultsOutImpuls$OOB_MAE

[1] 0.4119321

resultsOutImpuls$Rsquared

[1] 0.8901899

Train multiple seeds model

resultsOutImpuls <- multi_crfReg(dataIn=stimDataNum, iVars=iVars, dVar=dVar, seeds=seeds, ntree=ntree, mtry=mtry, permImpCondThres=permImpCondThres, nperm=nperm, minsplit=minsplit, minbucket=minbucket)

# print model prediction results
resultsOutImpuls$OOB_RMSE

[1] 0.531519

resultsOutImpuls$OOB_MAE

[1] 0.4075172

resultsOutImpuls$Rsquared

[1] 0.891422

# store results
resdAnnoyMnFitAB['Diff impuls', 'RMSE'] <- resultsOutImpuls$OOB_RMSE
resdAnnoyMnFitAB['Diff impuls', 'MAE'] <- resultsOutImpuls$OOB_MAE
resdAnnoyMnFitAB['Diff impuls', 'Rsquared'] <- resultsOutImpuls$Rsquared
resdAnnoyMnPermImpAB$DiffImpuls <- resultsOutImpuls$conditional_permimp

Plot results

par(mai=c(0,3,0,0))

# plot conditional importance
resultsOutImpuls.conimp <- arrange(resultsOutImpuls$conditional_permimp, desc(row_number()))

pBar <- ggplot(resultsOutImpuls.conimp) + geom_col(aes(x=factor(rownames(resultsOutImpuls.conimp), levels=rownames(resultsOutImpuls.conimp)), y=CondPermImp), fill=mycolours[6], width=0.5) + labs(x="Variable", y="Conditional variable permutation importance (mean change in annoyance)") + ggtitle("dImpulsiveness") + theme(text = element_text(family = "serif"), panel.grid=element_line(color = rgb(235, 235, 235, 100, maxColorValue = 255), linewidth = 0.25, linetype = 2)) + coord_flip()
pBar

if (saveplots){
  ggsave(filename="PtsABdAnnoyMndImpulsConPermimp.svg", width=8, height=3.6, path=file.path(outFigPath, "svg"))
  unlink("PtsABdAnnoyMndImpulsConPermimp.svg")
  
  ggsave(filename="PtsABdAnnoyMndImpulsConPermimp.pdf", width=8, height=3.6, path=file.path(outFigPath, "pdf"))
  unlink("PtsABdAnnoyMndImpulsConPermimp.pdf")
}

Selected metric

dImpulsVar <- "dImpulsLoudWZAvgMaxLR"

dSQM and loudness comparison

Now the highest importance dSQMs are ranked against each other, controlling for loudness difference.

Include dtonal loudness

Set variables

iVars <- c(diffVar, eventVar, dSharpVar, dTonLdVar, dFluctVar, dRoughVar, dImpulsVar)
dVar <- "dAnnoyMean"

seeds <- c(98465, 54163, 6541, 36485, 849675)

Hyperparameter tuning

p <- mtryTune(dataIn=stimDataNum, iVars=iVars, dVar=dVar, seeds=seeds[1:2],
             ntrees=ntrees, minsplit=minsplit, minbucket=minbucket)

p

Selected hyperparameters

ntree <- 501
mtry <- as.integer(length(iVars)/1.75)

Run model

Train preliminary model

nperm <- 5

resultsOutSQMs1 <- crfReg(dataIn=stimDataNum, iVars=iVars, dVar=dVar, seeds=seeds[1:2], ntree=ntree, mtry=mtry, permImpCondThres=permImpCondThres, nperm=nperm, minsplit=minsplit, minbucket=minbucket)

# print model prediction results
resultsOutSQMs1$OOB_RMSE

[1] 0.5266847

resultsOutSQMs1$OOB_MAE

[1] 0.4019546

resultsOutSQMs1$Rsquared

[1] 0.892996

Train multiple seeds model

resultsOutSQMs1 <- multi_crfReg(dataIn=stimDataNum, iVars=iVars, dVar=dVar, seeds=seeds, ntree=ntree, mtry=mtry, permImpCondThres=permImpCondThres, nperm=nperm, minsplit=minsplit, minbucket=minbucket)

# print model prediction results
resultsOutSQMs1$OOB_RMSE

[1] 0.5304116

resultsOutSQMs1$OOB_MAE

[1] 0.4046302

resultsOutSQMs1$Rsquared

[1] 0.8909152

# store results
resdAnnoyMnFitAB['Diff SQMs inc tonal loud', 'RMSE'] <- resultsOutSQMs1$OOB_RMSE
resdAnnoyMnFitAB['Diff SQMs inc tonal loud', 'MAE'] <- resultsOutSQMs1$OOB_MAE
resdAnnoyMnFitAB['Diff SQMs inc tonal loud', 'Rsquared'] <- resultsOutSQMs1$Rsquared
resdAnnoyMnPermImpAB$DiffSQMs1 <- resultsOutSQMs1$conditional_permimp

Plot results

par(mai=c(0,3,0,0))

# plot conditional importance
resultsOutSQMs1.conimp <- arrange(resultsOutSQMs1$conditional_permimp, desc(row_number()))

pBar <- ggplot(resultsOutSQMs1.conimp) + geom_col(aes(x=factor(rownames(resultsOutSQMs1.conimp), levels=rownames(resultsOutSQMs1.conimp)), y=CondPermImp), fill=mycolours[7], width=0.5) + labs(x="Variable", y="Conditional variable permutation importance (mean change in annoyance)") + theme(text = element_text(family = "serif"), panel.grid=element_line(color = rgb(235, 235, 235, 100, maxColorValue = 255), linewidth = 0.25, linetype = 2)) + coord_flip(ylim=c(0, 2.4))
pBar

if (saveplots){
  ggsave(filename="PtsABdAnnoyMnDiffSQMsTonLdConPermimp.svg", width=8, height=2.4, path=file.path(outFigPath, "svg"))
  unlink("PtsABdAnnoyMnDiffSQMsTonLdConPermimp.svg")
  
  ggsave(filename="PtsABdAnnoyMnDiffSQMsTonLdConPermimp.pdf", width=8, height=2.4, path=file.path(outFigPath, "pdf"))
  unlink("PtsABdAnnoyMnDiffSQMsTonLdConPermimp.pdf")
}

Exclude tonal loudness

Set variables

iVars <- c(diffVar, eventVar, dSharpVar, dTonalVar, dFluctVar, dRoughVar, dImpulsVar)
dVar <- "dAnnoyMean"

seeds <- c(49865, 7852, 845961, 410583, 36748)

Hyperparameter tuning

p <- mtryTune(dataIn=stimDataNum, iVars=iVars, dVar=dVar, seeds=seeds[1:2],
             ntrees=ntrees, minsplit=minsplit, minbucket=minbucket)

p

Selected hyperparameters

ntree <- 1501
mtry <- as.integer(length(iVars)/1.75)

Run model

Train preliminary model

nperm <- 5

resultsOutSQMs2 <- crfReg(dataIn=stimDataNum, iVars=iVars, dVar=dVar, seeds=seeds[1:2], ntree=ntree, mtry=mtry, permImpCondThres=permImpCondThres, nperm=nperm, minsplit=minsplit, minbucket=minbucket)

# print model prediction results
resultsOutSQMs2$OOB_RMSE

[1] 0.5229454

resultsOutSQMs2$OOB_MAE

[1] 0.4034912

resultsOutSQMs2$Rsquared

[1] 0.8960873

Train multiple seeds model

resultsOutSQMs2 <- multi_crfReg(dataIn=stimDataNum, iVars=iVars, dVar=dVar, seeds=seeds, ntree=ntree, mtry=mtry, permImpCondThres=permImpCondThres, nperm=nperm, minsplit=minsplit, minbucket=minbucket)

# print model prediction results
resultsOutSQMs2$OOB_RMSE

[1] 0.5230026

resultsOutSQMs2$OOB_MAE

[1] 0.4036246

resultsOutSQMs2$Rsquared

[1] 0.8962539

# store results
resdAnnoyMnFitAB['Diff SQMs no tonal loud', 'RMSE'] <- resultsOutSQMs2$OOB_RMSE
resdAnnoyMnFitAB['Diff SQMs no tonal loud', 'MAE'] <- resultsOutSQMs2$OOB_MAE
resdAnnoyMnFitAB['Diff SQMs no tonal loud', 'Rsquared'] <- resultsOutSQMs2$Rsquared
resdAnnoyMnPermImpAB$DiffSQMs2 <- resultsOutSQMs2$conditional_permimp

Plot results

par(mai=c(0,3,0,0))

# plot conditional importance
resultsOutSQMs2.conimp <- arrange(resultsOutSQMs2$conditional_permimp, desc(row_number()))

pBar <- ggplot(resultsOutSQMs2.conimp) + geom_col(aes(x=factor(rownames(resultsOutSQMs2.conimp), levels=rownames(resultsOutSQMs2.conimp)), y=CondPermImp), fill=mycolours[7], width=0.5) + labs(x="Variable", y="Conditional variable permutation importance (mean change in annoyance)") + theme(text = element_text(family = "serif"), panel.grid=element_line(color = rgb(235, 235, 235, 100, maxColorValue = 255), linewidth = 0.25, linetype = 2)) + coord_flip(ylim=c(0, 2.4))
pBar

if (saveplots){
  ggsave(filename="PtsABdAnnoyMnDiffSQMsNoTonLdConPermimp.svg", width=8, height=2.4, path=file.path(outFigPath, "svg"))
  unlink("PtsABdAnnoyMnDiffSQMsNoTonLdConPermimp.svg")
  
  ggsave(filename="PtsABdAnnoyMnDiffSQMsNoTonLdConPermimp.pdf", width=8, height=2.4, path=file.path(outFigPath, "pdf"))
  unlink("PtsABdAnnoyMnDiffSQMsNoTonLdConPermimp.pdf")
}

dPsychoacoustic annoyance metrics

Set variables

iVars <- c(ambVar, "dPsychAnnoyWidmann", "dPsychAnnoyMore", "dPsychAnnoyDi", "dPsychAnnoyTorija", "dPsychAnnoyWillemsen", "dPsychAnnoyBoucher")
dVar <- "dAnnoyMean"

seeds <- c(47896643, 475, 654, 98987132, 5446)

Hyperparameter tuning

p <- mtryTune(dataIn=stimDataNum, iVars=iVars, dVar=dVar, seeds=seeds[1:2],
              ntrees=ntrees, minsplit=minsplit, minbucket=minbucket)

p

Selected hyperparameters

ntree <- 1001
mtry <-  4

Run model

Train preliminary model

nperm <- 5

resultsOutPA <- crfReg(dataIn=stimDataNum, iVars=iVars, dVar=dVar, seeds=seeds[1:2], ntree=ntree, mtry=mtry, permImpCondThres=permImpCondThres, nperm=nperm, minsplit=minsplit, minbucket=minbucket)

# print model prediction results
resultsOutPA$OOB_RMSE

[1] 0.7781065

resultsOutPA$OOB_MAE

[1] 0.5943964

resultsOutPA$Rsquared

[1] 0.7605454

Train multiple seeds model

resultsOutPA <- multi_crfReg(dataIn=stimDataNum, iVars=iVars, dVar=dVar, seeds=seeds, ntree=ntree, mtry=mtry, permImpCondThres=permImpCondThres, nperm=nperm, minsplit=minsplit, minbucket=minbucket)

# print model prediction results
resultsOutPA$OOB_RMSE

[1] 0.7786348

resultsOutPA$OOB_MAE

[1] 0.5943763

resultsOutPA$Rsquared

[1] 0.7600659

# store results
resdAnnoyMnFitAB['Diff Psychoacoustic annoyance', 'RMSE'] <- resultsOutPA$OOB_RMSE
resdAnnoyMnFitAB['Diff Psychoacoustic annoyance', 'MAE'] <- resultsOutPA$OOB_MAE
resdAnnoyMnFitAB['Diff Psychoacoustic annoyance', 'Rsquared'] <- resultsOutPA$Rsquared
resdAnnoyMnPermImpAB$AbsPA <- resultsOutPA$conditional_permimp

Plot results

par(mai=c(0,3,0,0))

# plot conditional importance
resultsOutPA.conimp <- arrange(resultsOutPA$conditional_permimp, desc(row_number()))

pBar <- ggplot(resultsOutPA.conimp) + geom_col(aes(x=factor(rownames(resultsOutPA.conimp), levels=rownames(resultsOutPA.conimp)), y=CondPermImp), fill=mycolours[10], width=0.5) + labs(x="Variable", y="Conditional variable permutation importance (mean change in annoyance)") + theme(text = element_text(family = "serif"), panel.grid=element_line(color = rgb(235, 235, 235, 100, maxColorValue = 255), linewidth = 0.25, linetype = 2)) + coord_flip(ylim=c(0, 1.8))
pBar

if (saveplots){
  ggsave(filename="PtsABdAnnoyMndPAConPermimp.svg", width=8, height=2.4, path=file.path(outFigPath, "svg"))
  unlink("PtsABdAnnoyMndPAConPermimp.svg")
  
  ggsave(filename="PtsABdAnnoyMndPAConPermimp.pdf", width=8, height=2.4, path=file.path(outFigPath, "pdf"))
  unlink("PtsABdAnnoyMndPAConPermimp.pdf")
}

Save the results outputs to file

if (savedata){
  utils::write.csv(resdAnnoyMnFitAB, paste(outDataPath, "\\PtsABCRFdAnnoyMnOOBFit.csv", sep=""))
  ii <- 0
  temp = list()
  for (res in resdAnnoyMnPermImpAB){
    ii <- ii + 1
    temp[[ii]] <- as.data.frame(resdAnnoyMnPermImpAB[ii])
    names(temp[[ii]]) <- names(resdAnnoyMnPermImpAB[ii])
  }
  openxlsx::write.xlsx(temp, paste(outDataPath, "\\PtsABCRFdAnnoyMnConPermimp.xlsx",
                                   sep=""),
                       rowNames=TRUE)
}

(Change to) High annoyance

Initialise results output variables

resdHiAnnoyFitAB <- data.frame(RMSE = numeric(),
                             MAE = numeric(),
                             Rsquared = numeric())
resdHiAnnoyPermImpAB <- list()

All variables (absolute and difference)

Set variables

iVars <- names(stimDataNum)[which(names(stimDataNum) == 'UASEvents'):which(names(stimDataNum) == 'UASPsychAnnoyBoucher')]
iVars <- iVars[! iVars %in% 'SNRlevel']
iVars <- c(iVars,
           names(stimDataNum)[which(colnames(stimDataNum)=='LAeqLAF90diff'):
                               which(colnames(stimDataNum)=='dPsychAnnoyBoucher')], 'SNRlevel')
dVar <- "dHighAnnoyPc"

seeds <- c(2, 312, 1897, 465978, 821659)

Hyperparameter tuning

p <- mtryTune(dataIn=stimDataNum, iVars=iVars, dVar=dVar, seeds=seeds[1:2],
              ntrees=ntrees, minsplit=minsplit, minbucket=minbucket)

p

Selected hyperparameters

ntree <- 1501
mtry <- as.integer(length(iVars)/3.5)

Run model

Train preliminary model

nperm <- 5

resultsOutAbsDiffs <- crfReg(dataIn=stimDataNum, iVars=iVars, dVar=dVar, seeds=seeds[1:2], ntree=ntree, mtry=mtry, permImpCondThres=permImpCondThres, nperm=nperm, minsplit=minsplit, minbucket=minbucket)

# print model prediction results
resultsOutAbsDiffs$OOB_RMSE

[1] 6.145585

resultsOutAbsDiffs$OOB_MAE

[1] 4.687502

resultsOutAbsDiffs$Rsquared

[1] 0.6754223

Train multiple seeds model

resultsOutAbsDiffs <- multi_crfReg(dataIn=stimDataNum, iVars=iVars, dVar=dVar, seeds=seeds, ntree=ntree, mtry=mtry, permImpCondThres=permImpCondThres, nperm=nperm, minsplit=minsplit, minbucket=minbucket)

# print model prediction results
resultsOutAbsDiffs$OOB_RMSE

[1] 6.145665

resultsOutAbsDiffs$OOB_MAE

[1] 4.694338

resultsOutAbsDiffs$Rsquared

[1] 0.676116

# store results
resdHiAnnoyFitAB['All vars', 'RMSE'] <- resultsOutAbsDiffs$OOB_RMSE
resdHiAnnoyFitAB['All vars', 'MAE'] <- resultsOutAbsDiffs$OOB_MAE
resdHiAnnoyFitAB['All vars', 'Rsquared'] <- resultsOutAbsDiffs$Rsquared
resdHiAnnoyPermImpAB$AllVars <- resultsOutAbsDiffs$conditional_permimp

Plot results

par(mai=c(0,3,0,0))

# plot conditional importance
resultsOutAbsDiffs.conimp <- arrange(resultsOutAbsDiffs$conditional_permimp, desc(row_number()))

pBar <- ggplot(resultsOutAbsDiffs.conimp) + geom_col(aes(x=factor(rownames(resultsOutAbsDiffs.conimp), levels=rownames(resultsOutAbsDiffs.conimp)), y=CondPermImp), fill=mycolours[9], width=0.5) + labs(x="Variable", y="Conditional variable permutation importance (% highly annoyed)") + theme(text = element_text(family = "serif"), panel.grid=element_line(color = rgb(235, 235, 235, 100, maxColorValue = 255), linewidth = 0.25, linetype = 2)) + coord_flip()
pBar

if (saveplots){
  ggsave(filename="PtsABdHiAnnoyAllVarsConPermimp.svg", width=8, height=26, path=file.path(outFigPath, "svg"))
  unlink("PtsABdHiAnnoyAllVarsConPermimp.svg")

  ggsave(filename="PtsABdHiAnnoyAllVarsConPermimp.pdf", width=8, height=26, path=file.path(outFigPath, "pdf"))
  unlink("PtsABdHiAnnoyAllVarsConPermimp.pdf")
}

# Plot only positive values

resultsOutAbsDiffs.conimpPtv <- resultsOutAbsDiffs.conimp |>
                                          rownames_to_column('Metric') |>
                                                filter_if(is.numeric, all_vars(. > 0)) |>
                                                      column_to_rownames('Metric')

pBar <- ggplot(resultsOutAbsDiffs.conimpPtv) + geom_col(aes(x=factor(rownames(resultsOutAbsDiffs.conimpPtv), levels=rownames(resultsOutAbsDiffs.conimpPtv)), y=CondPermImp), fill=mycolours[9], width=0.5) + labs(x="Variable", y="Conditional variable permutation importance (% highly annoyed)") + theme(text = element_text(family = "serif"), panel.grid=element_line(color = rgb(235, 235, 235, 100, maxColorValue = 255), linewidth = 0.25, linetype = 2)) + coord_flip()
pBar

if (saveplots){
  ggsave(filename="PtsABdHiAnnoyAllVarsConPermimpPtv.svg", width=8, height=22, path=file.path(outFigPath, "svg"))
  unlink("PtsABdHiAnnoyAllVarsConPermimp.svg")
  
  ggsave(filename="PtsABdHiAnnoyAllVarsConPermimpPtv.pdf", width=8, height=22, path=file.path(outFigPath, "pdf"))
  unlink("PtsABdHiAnnoyAllVarsConPermimp.pdf")
}

# Plot only values within 1% of the maximum

resultsOutAbsDiffs.conimp1pc <- resultsOutAbsDiffs.conimp |>
                                          rownames_to_column('Metric') |>
                                                filter_if(is.numeric, all_vars(. > max(resultsOutAbsDiffs.conimp)/100)) |>
                                                      column_to_rownames('Metric')

pBar <- ggplot(resultsOutAbsDiffs.conimp1pc) + geom_col(aes(x=factor(rownames(resultsOutAbsDiffs.conimp1pc), levels=rownames(resultsOutAbsDiffs.conimp1pc)), y=CondPermImp), fill=mycolours[9], width=0.5) + labs(x="Variable", y="Conditional variable permutation importance (% highly annoyed)") + theme(text = element_text(family = "serif"), panel.grid=element_line(color = rgb(235, 235, 235, 100, maxColorValue = 255), linewidth = 0.25, linetype = 2)) + coord_flip()
pBar

if (saveplots){
  ggsave(filename="PtsABdHiAnnoyAllVarsConPermimp1pc.svg", width=8, height=7, path=file.path(outFigPath, "svg"))
  unlink("PtsABdHiAnnoyAllVarsConPermimp1pc.svg")
  
  ggsave(filename="PtsABdHiAnnoyAllVarsConPermimp1pc.pdf", width=8, height=7, path=file.path(outFigPath, "pdf"))
  unlink("PtsABdHiAnnoyAllVarsConPermimp1pc.pdf")
}

Absolute variables

Set variables

iVars <- names(stimDataNum)[which(names(stimDataNum) == 'UASEvents'):which(names(stimDataNum) == 'UASPsychAnnoyBoucher')]
iVars <- iVars[! iVars %in% c('SNRlevel', 'IntermitRatioC2MaxLR', 'IntermitRatioC3MaxLR', 'IntermitRatioC5MaxLR')]
dVar <- "dHighAnnoyPc"

seeds <- c(578312, 544, 84894, 54654, 153157)

Hyperparameter tuning

p <- mtryTune(dataIn=stimDataNum, iVars=iVars, dVar=dVar, seeds=seeds[1:2],
             ntrees=ntrees, minsplit=minsplit, minbucket=minbucket)

p

if (saveplots){
  ggsave(filename="PtsABdHiAnnoyAbsVarsHyperTune.svg", width=12, height=4, path=file.path(outFigPath, "svg"))
  unlink("PtsABdHiAnnoyAbsVarsHyperTune.svg")

  ggsave(filename="PtsABdHiAnnoyAbsVarsHyperTune.pdf", width=12, height=4, path=file.path(outFigPath, "pdf"))
  unlink("PtsABdHiAnnoyAbsVarsHyperTune.pdf")
}

Selected hyperparameters

ntree <- 1501
mtry <- as.integer(length(iVars)/1.75)

Run model

Train preliminary model

nperm <- 5

resultsOutAbs <- crfReg(dataIn=stimDataNum, iVars=iVars, dVar=dVar, seeds=seeds[1:2], ntree=ntree, mtry=mtry, permImpCondThres=permImpCondThres, nperm=nperm, minsplit=minsplit, minbucket=minbucket)

# print model prediction results
resultsOutAbs$OOB_RMSE

[1] 6.203493

resultsOutAbs$OOB_MAE

[1] 4.787191

resultsOutAbs$Rsquared

[1] 0.6664964

Train multiple seeds model

resultsOutAbs <- multi_crfReg(dataIn=stimDataNum, iVars=iVars, dVar=dVar, seeds=seeds, ntree=ntree, mtry=mtry, permImpCondThres=permImpCondThres, nperm=nperm, minsplit=minsplit, minbucket=minbucket)

# print model prediction results
resultsOutAbs$OOB_RMSE

[1] 6.205791

resultsOutAbs$OOB_MAE

[1] 4.775969

resultsOutAbs$Rsquared

[1] 0.6664051

# store results
resdHiAnnoyFitAB['Abs vars', 'RMSE'] <- resultsOutAbs$OOB_RMSE
resdHiAnnoyFitAB['Abs vars', 'MAE'] <- resultsOutAbs$OOB_MAE
resdHiAnnoyFitAB['Abs vars', 'Rsquared'] <- resultsOutAbs$Rsquared
resdHiAnnoyPermImpAB$AbsVars <- resultsOutAbs$conditional_permimp

Plot results

par(mai=c(0,3,0,0))

# plot conditional importance
resultsOutAbs.conimp <- arrange(resultsOutAbs$conditional_permimp, desc(row_number()))

pBar <- ggplot(resultsOutAbs.conimp) + geom_col(aes(x=factor(rownames(resultsOutAbs.conimp), levels=rownames(resultsOutAbs.conimp)), y=CondPermImp), fill=mycolours[1], width=0.5) + labs(x="Variable", y="Conditional variable permutation importance (% highly annoyed)") + theme(text = element_text(family = "serif"), panel.grid=element_line(color = rgb(235, 235, 235, 100, maxColorValue = 255), linewidth = 0.25, linetype = 2)) +
  coord_flip()
pBar

if (saveplots){
  ggsave(filename="PtsABdHiAnnoyAbsVarsConPermimp.svg", width=8, height=14, path=file.path(outFigPath, "svg"))
  unlink("PtsABdHiAnnoyAbsVarsConPermimp.svg")
  
  ggsave(filename="PtsABdHiAnnoyAbsVarsConPermimp.pdf", width=8, height=14, path=file.path(outFigPath, "pdf"))
  unlink("PtsABdHiAnnoyAbsVarsConPermimp.pdf")
}

# Plot only positive values
resultsOutAbs.conimpPtv <- resultsOutAbs.conimp |>
                                          rownames_to_column('Metric') |>
                                                filter_if(is.numeric, all_vars(. > 0)) |>
                                                      column_to_rownames('Metric')

pBar <- ggplot(resultsOutAbs.conimpPtv,) + geom_col(aes(x=factor(rownames(resultsOutAbs.conimpPtv), levels=rownames(resultsOutAbs.conimpPtv)), y=CondPermImp), fill=mycolours[1], width=0.5) + labs(x="Variable", y="Conditional variable permutation importance (% highly annoyed)") + theme(text = element_text(family = "serif"), panel.grid=element_line(color = rgb(235, 235, 235, 100, maxColorValue = 255), linewidth = 0.25, linetype = 2)) + coord_flip()
pBar

if (saveplots){
  ggsave(filename="PtsABdHiAnnoyAbsVarsConPermimpPtv.svg", width=8, height=10, path=file.path(outFigPath, "svg"))
  unlink("PtsABdHiAnnoyAbsVarsConPermimpPtv.svg")
  
  ggsave(filename="PtsABdHiAnnoyAbsVarsConPermimpPtv.pdf", width=8, height=10, path=file.path(outFigPath, "pdf"))
  unlink("PtsABdHiAnnoyAbsVarsConPermimpPtv.pdf")
}

# Plot only values within 1% of the maximum
resultsOutAbs.conimp1pc <- resultsOutAbs.conimp |>
                                          rownames_to_column('Metric') |>
                                                filter_if(is.numeric, all_vars(. > max(resultsOutAbs.conimp)/100)) |>
                                                      column_to_rownames('Metric')

pBar <- ggplot(resultsOutAbs.conimp1pc,) + geom_col(aes(x=factor(rownames(resultsOutAbs.conimp1pc), levels=rownames(resultsOutAbs.conimp1pc)), y=CondPermImp), fill=mycolours[1], width=0.5) + labs(x="Variable", y="Conditional variable permutation importance (% highly annoyed)") + theme(text = element_text(family = "serif"), panel.grid=element_line(color = rgb(235, 235, 235, 100, maxColorValue = 255), linewidth = 0.25, linetype = 2)) + coord_flip()
pBar

if (saveplots){
  ggsave(filename="PtsABdHiAnnoyAbsVarsConPermimp1pc.svg", width=8, height=4, path=file.path(outFigPath, "svg"))
  unlink("PtsABdHiAnnoyAbsVarsConPermimp1pc.svg")
  
  ggsave(filename="PtsABdHiAnnoyAbsVarsConPermimp1pc.pdf", width=8, height=4, path=file.path(outFigPath, "pdf"))
  unlink("PtsABdHiAnnoyAbsVarsConPermimp1pc.pdf")
}

Selected metric

absVar <- "UASLoudECMAPowAvgBin"

SQM analysis

Individual SQMs

Sharpness

Set variables

iVars <- c(absVar, eventVar, ambVar, "UASSharpAurISO3PowAvgBin", "UASSharpAurISO305ExBin", "UASSharpAurSHMPowAvgBin", "UASSharpAurSHM05ExBin", "UASSharpAurISO1PowAvgBin", "UASSharpAurISO105ExBin", "UASSharpvBISO1PowAvgBin", "UASSharpvBISO105ExBin", "UASSharpDINPowAvgBin", "UASSharpDIN05ExBin", "UASSharpAurISO1MedBin",
         "UASTonShpAurSHMPowAvgBin", "UASTonShpAurSHM05ExBin")
dVar <- "dHighAnnoyPc"

seeds <- c(7041, 905, 4984651, 6513213, 120651)

Hyperparameter tuning

p <- mtryTune(dataIn=stimDataNum, iVars=iVars, dVar=dVar, seeds=seeds[1:2],
             ntrees=ntrees, minsplit=minsplit, minbucket=minbucket)

p

NA
NA

Selected hyperparameters

ntree <- 1501
mtry <- as.integer(length(iVars)/2.25)

Run model

Train preliminary model

nperm <- 5

resultsOutSharp <- crfReg(dataIn=stimDataNum, iVars=iVars, dVar=dVar, seeds=seeds[1:2], ntree=ntree, mtry=mtry, permImpCondThres=permImpCondThres, nperm=nperm, minsplit=minsplit, minbucket=minbucket)

# print model prediction results
resultsOutSharp$OOB_RMSE

[1] 5.965448

resultsOutSharp$OOB_MAE

[1] 4.532257

resultsOutSharp$Rsquared

[1] 0.6957333

Train multiple seeds model

resultsOutSharp <- multi_crfReg(dataIn=stimDataNum, iVars=iVars, dVar=dVar, seeds=seeds, ntree=ntree, mtry=mtry, permImpCondThres=permImpCondThres, nperm=nperm, minsplit=minsplit, minbucket=minbucket)

# print model prediction results
resultsOutSharp$OOB_RMSE

[1] 5.945878

resultsOutSharp$OOB_MAE

[1] 4.51012

resultsOutSharp$Rsquared

[1] 0.6984185

# store results
resdHiAnnoyFitAB['Abs sharp', 'RMSE'] <- resultsOutSharp$OOB_RMSE
resdHiAnnoyFitAB['Abs sharp', 'MAE'] <- resultsOutSharp$OOB_MAE
resdHiAnnoyFitAB['Abs sharp', 'Rsquared'] <- resultsOutSharp$Rsquared
resdHiAnnoyPermImpAB$AbsSharp <- resultsOutSharp$conditional_permimp

Plot results

par(mai=c(0,3,0,0))

# plot conditional importance
resultsOutSharp.conimp <- arrange(resultsOutSharp$conditional_permimp, desc(row_number()))

pBar <- ggplot(resultsOutSharp.conimp) + geom_col(aes(x=factor(rownames(resultsOutSharp.conimp), levels=rownames(resultsOutSharp.conimp)), y=CondPermImp), fill=mycolours[2], width=0.5) + labs(x="Variable", y="Conditional variable permutation importance (% highly annoyed)") + ggtitle("Sharpness") + theme(text = element_text(family = "serif"), panel.grid=element_line(color = rgb(235, 235, 235, 100, maxColorValue = 255), linewidth = 0.25, linetype = 2)) + coord_flip()
pBar

if (saveplots){
  ggsave(filename="PtsABdHiAnnoySharpConPermimp.svg", width=8, height=4.9, path=file.path(outFigPath, "svg"))
  unlink("PtsABdHiAnnoySharpConPermimp.svg")
  
  ggsave(filename="PtsABdHiAnnoySharpConPermimp.pdf", width=8, height=4.9, path=file.path(outFigPath, "pdf"))
  unlink("PtsABdHiAnnoySharpConPermimp.pdf")
}

Selected metric

sharpVar <- "UASSharpAurISO3PowAvgBin"

Tonal loudness and tonality

Set variables

iVars <- c(absVar, eventVar, ambVar, "UASTonalECMAAvgMaxLR", "UASTonalSHMInt05ExMaxLR", "UASTonalSHMIntAvgMaxLR", "UASTonalECMA05ExMaxLR", "UASTonalAwSHMAvgMaxLR", "UASTonalAwSHM05ExMaxLR",   "UASTonalAwSHMIntAvgMaxLR",     "UASTonalAwSHMInt05ExMaxLR", "UASTonLdECMAPowAvgBin", "UASTonLdECMA05ExBin", "UASTonalAurAvgMaxLR", "UASTonalAur05ExMaxLR", "UASTonalAur10ExMaxLR",
         "UASTonShpAurSHMPowAvgBin", "UASTonShpAurSHM05ExBin")
dVar <- "dHighAnnoyPc"

seeds <- c(540, 104798, 456464, 87331, 94564)

Hyperparameter tuning

p <- mtryTune(dataIn=stimDataNum, iVars=iVars, dVar=dVar, seeds=seeds[1:2],
             ntrees=ntrees, minsplit=minsplit, minbucket=minbucket)

p

Selected hyperparameters

ntree <- 251
mtry <- as.integer(length(iVars)/1.5)

Run model

Train preliminary model

# Tonality with tonal loudness

nperm <- 5

resultsOutTonal1 <- crfReg(dataIn=stimDataNum, iVars=iVars, dVar=dVar, seeds=seeds[1:2], ntree=ntree, mtry=mtry, permImpCondThres=permImpCondThres, nperm=nperm, minsplit=minsplit, minbucket=minbucket)

# print model prediction results
resultsOutTonal1$OOB_RMSE

[1] 6.358715

resultsOutTonal1$OOB_MAE

[1] 4.808657

resultsOutTonal1$Rsquared

[1] 0.6484513

Train multiple seeds model

# Tonality with tonal loudness

resultsOutTonal1 <- multi_crfReg(dataIn=stimDataNum, iVars=iVars, dVar=dVar, seeds=seeds, ntree=ntree, mtry=mtry, permImpCondThres=permImpCondThres, nperm=nperm, minsplit=minsplit, minbucket=minbucket)

# print model prediction results
resultsOutTonal1$OOB_RMSE

[1] 6.320428

resultsOutTonal1$OOB_MAE

[1] 4.782126

resultsOutTonal1$Rsquared

[1] 0.6527379

# store results
resdHiAnnoyFitAB['Abs tonal inc loud', 'RMSE'] <- resultsOutTonal1$OOB_RMSE
resdHiAnnoyFitAB['Abs tonal inc loud', 'MAE'] <- resultsOutTonal1$OOB_MAE
resdHiAnnoyFitAB['Abs tonal inc loud', 'Rsquared'] <- resultsOutTonal1$Rsquared
resdHiAnnoyPermImpAB$AbsTonal1 <- resultsOutTonal1$conditional_permimp

Plot results

par(mai=c(0,3,0,0))

# plot conditional importance
resultsOutTonal1.conimp <- arrange(resultsOutTonal1$conditional_permimp, desc(row_number()))

pBar <- ggplot(resultsOutTonal1.conimp) + geom_col(aes(x=factor(rownames(resultsOutTonal1.conimp), levels=rownames(resultsOutTonal1.conimp)), y=CondPermImp), fill=mycolours[3], width=0.5) + labs(x="Variable", y="Conditional variable permutation importance (% highly annoyed)") + ggtitle("Tonality inc. tonal loudness") + theme(text = element_text(family = "serif"), panel.grid=element_line(color = rgb(235, 235, 235, 100, maxColorValue = 255), linewidth = 0.25, linetype = 2)) + coord_flip(ylim=c(0, 110))
pBar

if (saveplots){
  ggsave(filename="PtsABdHiAnnoyTonalLdConPermimp.svg", width=8, height=4.4, path=file.path(outFigPath, "svg"))
  unlink("PtsABdHiAnnoyTonalLdConPermimp.svg")
  
  ggsave(filename="PtsABdHiAnnoyTonalLdConPermimp.pdf", width=8, height=4.4, path=file.path(outFigPath, "pdf"))
  unlink("PtsABdHiAnnoyTonalLdConPermimp.pdf")
}

Selected metric

tonLdVar <- "UASTonLdECMAPowAvgBin"

Tonality without tonal loudness

Set variables

iVars <- c(absVar, eventVar, ambVar, "UASTonalECMAAvgMaxLR", "UASTonalSHMInt05ExMaxLR", "UASTonalSHMIntAvgMaxLR", "UASTonalECMA05ExMaxLR", "UASTonalAwSHMAvgMaxLR", "UASTonalAwSHM05ExMaxLR",   "UASTonalAwSHMIntAvgMaxLR", "UASTonalAwSHMInt05ExMaxLR",    "UASTonalAurAvgMaxLR", "UASTonalAur05ExMaxLR", "UASTonalAur10ExMaxLR")
dVar <- "dHighAnnoyPc"

seeds <- c(156089, 5860, 10528, 89541, 4685146)

Hyperparameter tuning

p <- mtryTune(dataIn=stimDataNum, iVars=iVars, dVar=dVar, seeds=seeds[1:2],
             ntrees=ntrees, minsplit=minsplit, minbucket=minbucket)

p

Selected hyperparameters

ntree <- 251
mtry <- as.integer(length(iVars)/1.25)

Run model

Train preliminary model

# Tonality

nperm <- 5

resultsOutTonal2 <- crfReg(dataIn=stimDataNum, iVars=iVars, dVar=dVar, seeds=seeds[1:2],
                           ntree=ntree, mtry=mtry, permImpCondThres=permImpCondThres, nperm=nperm,
                           minsplit=minsplit, minbucket=minbucket)

# print model prediction results
resultsOutTonal2$OOB_RMSE

[1] 6.434733

resultsOutTonal2$OOB_MAE

[1] 4.888113

resultsOutTonal2$Rsquared

[1] 0.640617

Train multiple seeds model

# Tonality

resultsOutTonal2 <- multi_crfReg(dataIn=stimDataNum, iVars=iVars, dVar=dVar, seeds=seeds, ntree=ntree, mtry=mtry, permImpCondThres=permImpCondThres, nperm=nperm, minsplit=minsplit, minbucket=minbucket)

# print model prediction results
resultsOutTonal2$OOB_RMSE

[1] 6.454118

resultsOutTonal2$OOB_MAE

[1] 4.898967

resultsOutTonal2$Rsquared

[1] 0.6386236

# store results
resdHiAnnoyFitAB['Abs tonal no loud', 'RMSE'] <- resultsOutTonal2$OOB_RMSE
resdHiAnnoyFitAB['Abs tonal no loud', 'MAE'] <- resultsOutTonal2$OOB_MAE
resdHiAnnoyFitAB['Abs tonal no loud', 'Rsquared'] <- resultsOutTonal2$Rsquared
resdHiAnnoyPermImpAB$AbsTonal2 <- resultsOutTonal2$conditional_permimp

Plot results

par(mai=c(0,3,0,0))

# plot conditional importance
resultsOutTonal2.conimp <- arrange(resultsOutTonal2$conditional_permimp, desc(row_number()))

pBar <- ggplot(resultsOutTonal2.conimp) + geom_col(aes(x=factor(rownames(resultsOutTonal2.conimp), levels=rownames(resultsOutTonal2.conimp)), y=CondPermImp), fill=mycolours[3], width=0.5) + labs(x="Variable", y="Conditional variable permutation importance (% highly annoyed)") + ggtitle("Tonality w/o tonal loudness") + theme(text = element_text(family = "serif"), panel.grid=element_line(color = rgb(235, 235, 235, 100, maxColorValue = 255), linewidth = 0.25, linetype = 2)) + coord_flip(ylim=c(0, 110))
pBar

if (saveplots){
  ggsave(filename="PtsABdHiAnnoyTonalConPermimp.svg", width=8, height=3.8, path=file.path(outFigPath, "svg"))
  unlink("PtsABdHiAnnoyTonalConPermimp.svg")
  
  ggsave(filename="PtsABdHiAnnoyTonalConPermimp.pdf", width=8, height=3.8, path=file.path(outFigPath, "pdf"))
  unlink("PtsABdHiAnnoyTonalConPermimp.pdf")
}

Selected metric

tonalVar <- "UASTonalAwSHMInt05ExMaxLR"

Fluctuation strength

Set variables

# Fluctuation strength
iVars <- c(absVar, eventVar, ambVar, "UASFluctOldSHM10ExBin", "UASFluctOldSHM05ExBin", "UASFluctECMA10ExBin", "UASFluctECMA05ExBin", "UASFluctFZ10ExMaxLR", "UASFluctFZ05ExMaxLR", "UASFluctOV10ExMaxLR", "UASFluctOV05ExMaxLR")
dVar <- "dHighAnnoyPc"

seeds <- c(25107, 546098, 195, 5937, 102658)

Hyperparameter tuning

p <- mtryTune(dataIn=stimDataNum, iVars=iVars, dVar=dVar, seeds=seeds[1:2],
             ntrees=ntrees, minsplit=minsplit, minbucket=minbucket)

p

Selected hyperparameters

ntree <- 251
mtry <- as.integer(length(iVars)/1.25)

Run model

Train preliminary model

nperm <- 5

resultsOutFluct <- crfReg(dataIn=stimDataNum, iVars=iVars, dVar=dVar, seeds=seeds[1:2],
                          ntree=ntree, mtry=mtry, permImpCondThres=permImpCondThres,
                          nperm=nperm, minsplit=minsplit, minbucket=minbucket)

# print model prediction results
resultsOutFluct$OOB_RMSE

[1] 6.468882

resultsOutFluct$OOB_MAE

[1] 4.835973

resultsOutFluct$Rsquared

[1] 0.636221

Train multiple seeds model

resultsOutFluct <- multi_crfReg(dataIn=stimDataNum, iVars=iVars, dVar=dVar, seeds=seeds, ntree=ntree, mtry=mtry, permImpCondThres=permImpCondThres, nperm=nperm, minsplit=minsplit, minbucket=minbucket)

# print model prediction results
resultsOutFluct$OOB_RMSE

[1] 6.48319

resultsOutFluct$OOB_MAE

[1] 4.846364

resultsOutFluct$Rsquared

[1] 0.6345453

# store results
resdHiAnnoyFitAB['Abs fluct', 'RMSE'] <- resultsOutFluct$OOB_RMSE
resdHiAnnoyFitAB['Abs fluct', 'MAE'] <- resultsOutFluct$OOB_MAE
resdHiAnnoyFitAB['Abs fluct', 'Rsquared'] <- resultsOutFluct$Rsquared
resdHiAnnoyPermImpAB$AbsFluct <- resultsOutFluct$conditional_permimp

Plot results

par(mai=c(0,3,0,0))

# plot conditional importance
resultsOutFluct.conimp <- arrange(resultsOutFluct$conditional_permimp, desc(row_number()))

pBar <- ggplot(resultsOutFluct.conimp) + geom_col(aes(x=factor(rownames(resultsOutFluct.conimp), levels=rownames(resultsOutFluct.conimp)), y=CondPermImp), fill=mycolours[4], width=0.5) + labs(x="Variable", y="Conditional variable permutation importance (% highly annoyed)") + ggtitle("Fluctuation strength") + theme(text = element_text(family = "serif"), panel.grid=element_line(color = rgb(235, 235, 235, 100, maxColorValue = 255), linewidth = 0.25, linetype = 2)) + coord_flip()
pBar

if (saveplots){
  ggsave(filename="PtsABdHiAnnoyFluctConPermimp.svg", width=8, height=2.9, path=file.path(outFigPath, "svg"))
  unlink("PtsABdHiAnnoyFluctConPermimp.svg")
  
  ggsave(filename="PtsABdHiAnnoyFluctConPermimp.pdf", width=8, height=2.9, path=file.path(outFigPath, "pdf"))
  unlink("PtsABdHiAnnoyFluctConPermimp.pdf")
}

Selected metric

fluctVar <- "UASFluctECMA10ExBin"

Roughness

Set variables

# Roughness
iVars <- c(absVar, eventVar, ambVar, "UASRoughECMA10ExBin", "UASRoughECMA05ExBin", "UASRoughFZ10ExMaxLR", "UASRoughFZ05ExMaxLR", "UASRoughDW10ExMaxLR", "UASRoughDW05ExMaxLR")
dVar <- "dHighAnnoyPc"

seeds <- c(4701, 52187, 16589, 65217, 16893)

Hyperparameter tuning

p <- mtryTune(dataIn=stimDataNum, iVars=iVars, dVar=dVar, seeds=seeds[1:2],
             ntrees=ntrees, minsplit=minsplit, minbucket=minbucket)

p

Selected hyperparameters

ntree <- 4001
mtry <- as.integer(length(iVars)/1.5)

Run model

Train preliminary model

nperm <- 5

resultsOutRough <- crfReg(dataIn=stimDataNum, iVars=iVars, dVar=dVar, seeds=seeds[1:2], ntree=ntree, mtry=mtry, permImpCondThres=permImpCondThres, nperm=nperm, minsplit=minsplit, minbucket=minbucket)

# print model prediction results
resultsOutRough$OOB_RMSE

[1] 6.370114

resultsOutRough$OOB_MAE

[1] 4.823308

resultsOutRough$Rsquared

[1] 0.6515551

Train multiple seeds model

resultsOutRough <- multi_crfReg(dataIn=stimDataNum, iVars=iVars, dVar=dVar, seeds=seeds, ntree=ntree, mtry=mtry, permImpCondThres=permImpCondThres, nperm=nperm, minsplit=minsplit, minbucket=minbucket)

# print model prediction results
resultsOutRough$OOB_RMSE

[1] 6.398763

resultsOutRough$OOB_MAE

[1] 4.841106

resultsOutRough$Rsquared

[1] 0.6481992

# store results
resdHiAnnoyFitAB['Abs rough', 'RMSE'] <- resultsOutRough$OOB_RMSE
resdHiAnnoyFitAB['Abs rough', 'MAE'] <- resultsOutRough$OOB_MAE
resdHiAnnoyFitAB['Abs rough', 'Rsquared'] <- resultsOutRough$Rsquared
resdHiAnnoyPermImpAB$AbsRough <- resultsOutRough$conditional_permimp

Plot results

par(mai=c(0,3,0,0))

# plot conditional importance
resultsOutRough.conimp <- arrange(resultsOutRough$conditional_permimp, desc(row_number()))

pBar <- ggplot(resultsOutRough.conimp) + geom_col(aes(x=factor(rownames(resultsOutRough.conimp), levels=rownames(resultsOutRough.conimp)), y=CondPermImp), fill=mycolours[5], width=0.5) + labs(x="Variable", y="Conditional variable permutation importance (% highly annoyed)") + ggtitle("Roughness") + theme(text = element_text(family = "serif"), panel.grid=element_line(color = rgb(235, 235, 235, 100, maxColorValue = 255), linewidth = 0.25, linetype = 2)) + coord_flip()
pBar

if (saveplots){
  ggsave(filename="PtsABdHiAnnoyRoughConPermimp.svg", width=8, height=2.9, path=file.path(outFigPath, "svg"))
  unlink("PtsABdHiAnnoyRoughConPermimp.svg")
  
  ggsave(filename="PtsABdHiAnnoyRoughConPermimp.pdf", width=8, height=2.9, path=file.path(outFigPath, "pdf"))
  unlink("PtsABdHiAnnoyRoughConPermimp.pdf")
}

Selected metric

roughVar <- "UASRoughFZ05ExMaxLR"

Impulsiveness

Set variables

# Impulsiveness
iVars <- c(absVar, eventVar, ambVar, "UASImpulsSHMAvgMaxLR", "UASImpulsSHM05ExMaxLR", "UASImpulsSHMPowAvgMaxLR", "UASImpulsLoudWZAvgMaxLR", "UASImpulsLoudWZ05ExMaxLR", "UASImpulsLoudWZPowAvgMaxLR", "UASImpulsLoudWECMAAvgBin", "UASImpulsLoudWECMA05ExBin", "UASImpulsLoudWECMAPowAvgBin")
dVar <- "dHighAnnoyPc"

seeds <- c(8495, 59867, 5416, 9843, 86)

Hyperparameter tuning

p <- mtryTune(dataIn=stimDataNum, iVars=iVars, dVar=dVar, seeds=seeds[1:2],
             ntrees=ntrees, minsplit=minsplit, minbucket=minbucket)

p

Selected hyperparameters

ntree <- 501
mtry <- as.integer(length(iVars)/1.5)

Run model

Train preliminary model

nperm <- 5

resultsOutImpuls <- crfReg(dataIn=stimDataNum, iVars=iVars, dVar=dVar, seeds=seeds[1:2], ntree=ntree, mtry=mtry, permImpCondThres=permImpCondThres, nperm=nperm, minsplit=minsplit, minbucket=minbucket)

# print model prediction results
resultsOutImpuls$OOB_RMSE

[1] 6.270202

resultsOutImpuls$OOB_MAE

[1] 4.842259

resultsOutImpuls$Rsquared

[1] 0.6582782

Train multiple seeds model

resultsOutImpuls <- multi_crfReg(dataIn=stimDataNum, iVars=iVars, dVar=dVar, seeds=seeds, ntree=ntree, mtry=mtry, permImpCondThres=permImpCondThres, nperm=nperm, minsplit=minsplit, minbucket=minbucket)

# print model prediction results
resultsOutImpuls$OOB_RMSE

[1] 6.270233

resultsOutImpuls$OOB_MAE

[1] 4.84118

resultsOutImpuls$Rsquared

[1] 0.6582749

# store results
resdHiAnnoyFitAB['Abs impuls', 'RMSE'] <- resultsOutImpuls$OOB_RMSE
resdHiAnnoyFitAB['Abs impuls', 'MAE'] <- resultsOutImpuls$OOB_MAE
resdHiAnnoyFitAB['Abs impuls', 'Rsquared'] <- resultsOutImpuls$Rsquared
resdHiAnnoyPermImpAB$AbsImpuls <- resultsOutImpuls$conditional_permimp

Plot results

par(mai=c(0,3,0,0))

# plot conditional importance
resultsOutImpuls.conimp <- arrange(resultsOutImpuls$conditional_permimp, desc(row_number()))

pBar <- ggplot(resultsOutImpuls.conimp) + geom_col(aes(x=factor(rownames(resultsOutImpuls.conimp), levels=rownames(resultsOutImpuls.conimp)), y=CondPermImp), fill=mycolours[6], width=0.5) + labs(x="Variable", y="Conditional variable permutation importance (% highly annoyed)") + ggtitle("Impulsiveness") + theme(text = element_text(family = "serif"), panel.grid=element_line(color = rgb(235, 235, 235, 100, maxColorValue = 255), linewidth = 0.25, linetype = 2)) + coord_flip()
pBar

if (saveplots){
  ggsave(filename="PtsABdHiAnnoyImpulsConPermimp.svg", width=8, height=4, path=file.path(outFigPath, "svg"))
  unlink("PtsABdHiAnnoyImpulsConPermimp.svg")
  
  ggsave(filename="PtsABdHiAnnoyImpulsConPermimp.pdf", width=8, height=4, path=file.path(outFigPath, "pdf"))
  unlink("PtsABdHiAnnoyImpulsConPermimp.pdf")
}

Selected metric

impulsVar <- "UASImpulsLoudWZAvgMaxLR"

SQM and loudness comparison

Now the highest importance SQMs are ranked against each other, controlling for UAS loudness and ambient LAeq.

Include tonal loudness

Set variables

iVars <- c(absVar, eventVar, ambVar, sharpVar, tonLdVar, fluctVar, roughVar, impulsVar)
dVar <- "dHighAnnoyPc"

seeds <- c(70498, 4, 14986, 453, 864)

Hyperparameter tuning

p <- mtryTune(dataIn=stimDataNum, iVars=iVars, dVar=dVar, seeds=seeds[1:2],
             ntrees=ntrees, minsplit=minsplit, minbucket=minbucket)

p

Selected hyperparameters

ntree <- 1501
mtry <- 3

Run model

Train preliminary model

nperm <- 5

resultsOutSQMs1 <- crfReg(dataIn=stimDataNum, iVars=iVars, dVar=dVar, seeds=seeds[1:2], ntree=ntree, mtry=mtry, permImpCondThres=permImpCondThres, nperm=nperm, minsplit=minsplit, minbucket=minbucket)

# print model prediction results
resultsOutSQMs1$OOB_RMSE

[1] 6.112539

resultsOutSQMs1$OOB_MAE

[1] 4.697572

resultsOutSQMs1$Rsquared

[1] 0.6758278

Train multiple seeds model

resultsOutSQMs1 <- multi_crfReg(dataIn=stimDataNum, iVars=iVars, dVar=dVar, seeds=seeds, ntree=ntree, mtry=mtry, permImpCondThres=permImpCondThres, nperm=nperm, minsplit=minsplit, minbucket=minbucket)

# print model prediction results
resultsOutSQMs1$OOB_RMSE

[1] 6.115048

resultsOutSQMs1$OOB_MAE

[1] 4.695767

resultsOutSQMs1$Rsquared

[1] 0.6754898

# store results
resdHiAnnoyFitAB['Abs SQMs inc tonal loud', 'RMSE'] <- resultsOutSQMs1$OOB_RMSE
resdHiAnnoyFitAB['Abs SQMs inc tonal loud', 'MAE'] <- resultsOutSQMs1$OOB_MAE
resdHiAnnoyFitAB['Abs SQMs inc tonal loud', 'Rsquared'] <- resultsOutSQMs1$Rsquared
resdHiAnnoyPermImpAB$AbsSQMs1 <- resultsOutSQMs1$conditional_permimp

Plot results

par(mai=c(0,3,0,0))

# plot conditional importance
resultsOutSQMs1.conimp <- arrange(resultsOutSQMs1$conditional_permimp, desc(row_number()))

pBar <- ggplot(resultsOutSQMs1.conimp) + geom_col(aes(x=factor(rownames(resultsOutSQMs1.conimp), levels=rownames(resultsOutSQMs1.conimp)), y=CondPermImp), fill=mycolours[7], width=0.5) + labs(x="Variable", y="Conditional variable permutation importance (% highly annoyed)") + theme(text = element_text(family = "serif"), panel.grid=element_line(color = rgb(235, 235, 235, 100, maxColorValue = 255), linewidth = 0.25, linetype = 2)) + coord_flip(ylim=c(0, 90))
pBar

if (saveplots){
  ggsave(filename="PtsABdHiAnnoyAbsSQMsTonLdConPermimp.svg", width=8, height=2.4, path=file.path(outFigPath, "svg"))
  unlink("PtsABdHiAnnoyAbsSQMsTonLdConPermimp.svg")
  
  ggsave(filename="PtsABdHiAnnoyAbsSQMsTonLdConPermimp.pdf", width=8, height=2.4, path=file.path(outFigPath, "pdf"))
  unlink("PtsABdHiAnnoyAbsSQMsTonLdConPermimp.pdf")
}

Exclude tonal loudness

Set variables

iVars <- c(absVar, eventVar, ambVar, sharpVar, tonalVar, fluctVar, roughVar, impulsVar)
dVar <- "dHighAnnoyPc"

seeds <- c(546, 57203, 270835, 60592, 8094)

Hyperparameter tuning

p <- mtryTune(dataIn=stimDataNum, iVars=iVars, dVar=dVar, seeds=seeds[1:2],
             ntrees=ntrees, minsplit=minsplit, minbucket=minbucket)

p

Selected hyperparameters

ntree <- 4001
mtry <- 3

Run model

Train preliminary model

nperm <- 5

resultsOutSQMs2 <- crfReg(dataIn=stimDataNum, iVars=iVars, dVar=dVar, seeds=seeds[1:2], ntree=ntree, mtry=mtry, permImpCondThres=permImpCondThres, nperm=nperm, minsplit=minsplit, minbucket=minbucket)

# print model prediction results
resultsOutSQMs2$OOB_RMSE

[1] 6.081512

resultsOutSQMs2$OOB_MAE

[1] 4.686507

resultsOutSQMs2$Rsquared

[1] 0.6810178

Train multiple seeds model

resultsOutSQMs2 <- multi_crfReg(dataIn=stimDataNum, iVars=iVars, dVar=dVar, seeds=seeds, ntree=ntree, mtry=mtry, permImpCondThres=permImpCondThres, nperm=nperm, minsplit=minsplit, minbucket=minbucket)

# print model prediction results
resultsOutSQMs2$OOB_RMSE

[1] 6.101152

resultsOutSQMs2$OOB_MAE

[1] 4.702343

resultsOutSQMs2$Rsquared

[1] 0.6787822

# store results
resdHiAnnoyFitAB['Abs SQMs no tonal loud', 'RMSE'] <- resultsOutSQMs2$OOB_RMSE
resdHiAnnoyFitAB['Abs SQMs no tonal loud', 'MAE'] <- resultsOutSQMs2$OOB_MAE
resdHiAnnoyFitAB['Abs SQMs no tonal loud', 'Rsquared'] <- resultsOutSQMs2$Rsquared
resdHiAnnoyPermImpAB$AbsSQMs2 <- resultsOutSQMs2$conditional_permimp

Plot results

par(mai=c(0,3,0,0))

# plot conditional importance
resultsOutSQMs2.conimp <- arrange(resultsOutSQMs2$conditional_permimp, desc(row_number()))

pBar <- ggplot(resultsOutSQMs2.conimp) + geom_col(aes(x=factor(rownames(resultsOutSQMs2.conimp), levels=rownames(resultsOutSQMs2.conimp)), y=CondPermImp), fill=mycolours[7], width=0.5) + labs(x="Variable", y="Conditional variable permutation importance (% highly annoyed)") + theme(text = element_text(family = "serif"), panel.grid=element_line(color = rgb(235, 235, 235, 100, maxColorValue = 255), linewidth = 0.25, linetype = 2)) + coord_flip(ylim=c(0, 90))
pBar

if (saveplots){
  ggsave(filename="PtsABdHiAnnoyAbsSQMsNoTonLdConPermimp.svg", width=8, height=2.4, path=file.path(outFigPath, "svg"))
  unlink("PtsABdHiAnnoyAbsSQMsNoTonLdConPermimp.svg")
  
  ggsave(filename="PtsABdHiAnnoyAbsSQMsNoTonLdConPermimp.pdf", width=8, height=2.4, path=file.path(outFigPath, "pdf"))
  unlink("PtsABdHiAnnoyAbsSQMsNoTonLdConPermimp.pdf")
}

Psychoacoustic annoyance metrics

Set variables

iVars <- c(ambVar, "UASPsychAnnoyWidmann", "UASPsychAnnoyMore", "UASPsychAnnoyDi", "UASPsychAnnoyTorija", "UASPsychAnnoyWillemsen", "UASPsychAnnoyBoucher")
dVar <- "dHighAnnoyPc"

seeds <- c(48651, 45, 785123, 65, 5163)

Hyperparameter tuning

p <- mtryTune(dataIn=stimDataNum, iVars=iVars, dVar=dVar, seeds=seeds[1:2],
             ntrees=ntrees, minsplit=minsplit, minbucket=minbucket)

p

Selected hyperparameters

ntree <- 4001
mtry <- 4

Run model

Train preliminary model

nperm <- 5

resultsOutPA <- crfReg(dataIn=stimDataNum, iVars=iVars, dVar=dVar, seeds=seeds[1:2], ntree=ntree, mtry=mtry, permImpCondThres=permImpCondThres, nperm=nperm, minsplit=minsplit, minbucket=minbucket)

# print model prediction results
resultsOutPA$OOB_RMSE

[1] 6.88072

resultsOutPA$OOB_MAE

[1] 5.141535

resultsOutPA$Rsquared

[1] 0.5875095

Train multiple seeds model

resultsOutPA <- multi_crfReg(dataIn=stimDataNum, iVars=iVars, dVar=dVar, seeds=seeds, ntree=ntree, mtry=mtry, permImpCondThres=permImpCondThres, nperm=nperm, minsplit=minsplit, minbucket=minbucket)

# print model prediction results
resultsOutPA$OOB_RMSE

[1] 6.879401

resultsOutPA$OOB_MAE

[1] 5.138723

resultsOutPA$Rsquared

[1] 0.5876636

# store results
resdHiAnnoyFitAB['Psychoacoustic annoyance', 'RMSE'] <- resultsOutPA$OOB_RMSE
resdHiAnnoyFitAB['Psychoacoustic annoyance', 'MAE'] <- resultsOutPA$OOB_MAE
resdHiAnnoyFitAB['Psychoacoustic annoyance', 'Rsquared'] <- resultsOutPA$Rsquared
resdHiAnnoyPermImpAB$AbsPA <- resultsOutPA$conditional_permimp

Plot results

par(mai=c(0,3,0,0))

# plot conditional importance
resultsOutPA.conimp <- arrange(resultsOutPA$conditional_permimp, desc(row_number()))

pBar <- ggplot(resultsOutPA.conimp) + geom_col(aes(x=factor(rownames(resultsOutPA.conimp), levels=rownames(resultsOutPA.conimp)), y=CondPermImp), fill=mycolours[10], width=0.5) + labs(x="Variable", y="Conditional variable permutation importance (% highly annoyed)") + theme(text = element_text(family = "serif"), panel.grid=element_line(color = rgb(235, 235, 235, 100, maxColorValue = 255), linewidth = 0.25, linetype = 2)) + coord_flip(ylim=c(0, 60))
pBar

if (saveplots){
  ggsave(filename="PtsABdHiAnnoyAbsPAConPermimp.svg", width=8, height=2.4, path=file.path(outFigPath, "svg"))
  unlink("PtsABdHiAnnoyAbsPAConPermimp.svg")
  
  ggsave(filename="PtsABdHiAnnoyAbsPAConPermimp.pdf", width=8, height=2.4, path=file.path(outFigPath, "pdf"))
  unlink("PtsABdHiAnnoyAbsPAConPermimp.pdf")
}

Difference variables

Next, the difference metrics are analysed

Set variables

iVars <- c(names(stimDataNum)[which(colnames(stimDataNum)=="LAeqLAF90diff"):
                               which(colnames(stimDataNum)=="dPsychAnnoyBoucher")],
           "UASEvents", "SNRlevel", "IntermitRatioC2MaxLR", "IntermitRatioC3MaxLR", "IntermitRatioC5MaxLR")
dVar <- "dHighAnnoyPc"

seeds <- c(568392, 498, 4089, 78132, 741809)

Hyperparameter tuning

p <- mtryTune(dataIn=stimDataNum, iVars=iVars, dVar=dVar, seeds=seeds[1:2],
             ntrees=ntrees, minsplit=minsplit, minbucket=minbucket)

p

Selected hyperparameters

ntree <- 1501
mtry <- as.integer(length(iVars)/5.25)

Run model

Train preliminary model

nperm <- 5

resultsOutDiffs <- crfReg(dataIn=stimDataNum, iVars=iVars, dVar=dVar, seeds=seeds[1:2], ntree=ntree, mtry=mtry, permImpCondThres=permImpCondThres, nperm=nperm, minsplit=minsplit, minbucket=minbucket)

# print model prediction results
resultsOutDiffs$OOB_RMSE

[1] 6.924287

resultsOutDiffs$OOB_MAE

[1] 5.437204

resultsOutDiffs$Rsquared

[1] 0.5889082

Train multiple seeds model

resultsOutDiffs <- multi_crfReg(dataIn=stimDataNum, iVars=iVars, dVar=dVar, seeds=seeds, ntree=ntree, mtry=mtry, permImpCondThres=permImpCondThres, nperm=nperm, minsplit=minsplit, minbucket=minbucket)

# print model prediction results
resultsOutDiffs$OOB_RMSE

[1] 6.890232

resultsOutDiffs$OOB_MAE

[1] 5.40236

resultsOutDiffs$Rsquared

[1] 0.5926204

# store results
resdHiAnnoyFitAB['Diff vars', 'RMSE'] <- resultsOutDiffs$OOB_RMSE
resdHiAnnoyFitAB['Diff vars', 'MAE'] <- resultsOutDiffs$OOB_MAE
resdHiAnnoyFitAB['Diff vars', 'Rsquared'] <- resultsOutDiffs$Rsquared
resdHiAnnoyPermImpAB$DiffVars <- resultsOutDiffs$conditional_permimp

Plot results

par(mai=c(0,3,0,0))

# plot conditional importance
resultsOutDiffs.conimp <- arrange(resultsOutDiffs$conditional_permimp, desc(row_number()))

pBar <- ggplot(resultsOutDiffs.conimp) + geom_col(aes(x=factor(rownames(resultsOutDiffs.conimp), levels=rownames(resultsOutDiffs.conimp)), y=CondPermImp), fill=mycolours[8], width=0.5) + labs(x="Variable", y="Conditional variable permutation importance (% highly annoyed)") + theme(text = element_text(family = "serif"), panel.grid=element_line(color = rgb(235, 235, 235, 100, maxColorValue = 255), linewidth = 0.25, linetype = 2)) + coord_flip()
pBar

if (saveplots){
  ggsave(filename="PtsABdHiAnnoyDiffVarsConPermimp.svg", width=8, height=16, path=file.path(outFigPath, "svg"))
  unlink("PtsABdHiAnnoyDiffVarsConPermimp.svg")
  
  ggsave(filename="PtsABdHiAnnoyDiffVarsConPermimp.pdf", width=8, height=16, path=file.path(outFigPath, "pdf"))
  unlink("PtsABdHiAnnoyDiffVarsConPermimp.pdf")
}

# Plot only positive values
resultsOutDiffs.conimpPtv <- resultsOutDiffs.conimp |>
                                        rownames_to_column('Metric') |>
                                            filter_if(is.numeric, all_vars(. > 0)) |>
                                                  column_to_rownames('Metric')

pBar <- ggplot(resultsOutDiffs.conimpPtv) + geom_col(aes(x=factor(rownames(resultsOutDiffs.conimpPtv), levels=rownames(resultsOutDiffs.conimpPtv)), y=CondPermImp), fill=mycolours[8], width=0.5) + labs(x="Variable", y="Conditional variable permutation importance (% highly annoyed)") + theme(text = element_text(family = "serif"), panel.grid=element_line(color = rgb(235, 235, 235, 100, maxColorValue = 255), linewidth = 0.25, linetype = 2)) + coord_flip()
pBar

if (saveplots){
  ggsave(filename="PtsABdHiAnnoyDiffVarsConPermimpPtv.svg", width=8, height=14, path=file.path(outFigPath, "svg"))
  unlink("PtsABdHiAnnoyDiffVarsConPermimpPtv.svg")
  
  ggsave(filename="PtsABdHiAnnoyDiffVarsConPermimpPtv.pdf", width=8, height=14, path=file.path(outFigPath, "pdf"))
  unlink("PtsABdHiAnnoyDiffVarsConPermimpPtv.pdf")
}

# Plot only values within 1% of the maximum
resultsOutDiffs.conimp1pc <- resultsOutDiffs.conimp |>
                                        rownames_to_column('Metric') |>
                                            filter_if(is.numeric, all_vars(. > max(resultsOutDiffs.conimp)/100)) |>
                                                  column_to_rownames('Metric')

pBar <- ggplot(resultsOutDiffs.conimp1pc) + geom_col(aes(x=factor(rownames(resultsOutDiffs.conimp1pc), levels=rownames(resultsOutDiffs.conimp1pc)), y=CondPermImp), fill=mycolours[8], width=0.5) + labs(x="Variable", y="Conditional variable permutation importance (% highly annoyed)") + theme(text = element_text(family = "serif"), panel.grid=element_line(color = rgb(235, 235, 235, 100, maxColorValue = 255), linewidth = 0.25, linetype = 2)) + coord_flip()
pBar

if (saveplots){
  ggsave(filename="PtsABdHiAnnoyDiffVarsConPermimp1pc.svg", width=8, height=8, path=file.path(outFigPath, "svg"))
  unlink("PtsABdHiAnnoyDiffVarsConPermimp1pc.svg")
  
  ggsave(filename="PtsABdHiAnnoyDiffVarsConPermimp1pc.pdf", width=8, height=8, path=file.path(outFigPath, "pdf"))
  unlink("PtsABdHiAnnoyDiffVarsConPermimp1pc.pdf")
}

Selected metric

diffVar <- "EPNLLAeqdiff"

dSQM analysis

Individual SQMs

dSharpness

Set variables

iVars <- c(diffVar, eventVar, "dSharpAurISO3PowAvgBin", "dSharpAurISO305ExBin", "dSharpAurSHMPowAvgBin", "dSharpAurSHM05ExBin", "dTonShpAurSHMPowAvgBin", "dTonShpAurSHM05ExBin", "PartTonShpAurSHMPowAvgBin",
           "PartTonShpAurSHM05ExBin")
dVar <- "dHighAnnoyPc"

seeds <- c(84194, 905, 64815, 928054, 625091, 582031)

Hyperparameter tuning

p <- mtryTune(dataIn=stimDataNum, iVars=iVars, dVar=dVar, seeds=seeds[1:2],
             ntrees=ntrees, minsplit=minsplit, minbucket=minbucket)

p

Selected hyperparameters

ntree <- 251
mtry <- as.integer(length(iVars)/2.25)

Run model

Train preliminary model

nperm <- 5

resultsOutSharp <- crfReg(dataIn=stimDataNum, iVars=iVars, dVar=dVar, seeds=seeds[1:2], ntree=ntree, mtry=mtry, permImpCondThres=permImpCondThres, nperm=nperm, minsplit=minsplit, minbucket=minbucket)

# print model prediction results
resultsOutSharp$OOB_RMSE

[1] 6.684201

resultsOutSharp$OOB_MAE

[1] 5.128462

resultsOutSharp$Rsquared

[1] 0.6129086

Train multiple seeds model

resultsOutSharp <- multi_crfReg(dataIn=stimDataNum, iVars=iVars, dVar=dVar, seeds=seeds, ntree=ntree, mtry=mtry, permImpCondThres=permImpCondThres, nperm=nperm, minsplit=minsplit, minbucket=minbucket)

# print model prediction results
resultsOutSharp$OOB_RMSE

[1] 6.636169

resultsOutSharp$OOB_MAE

[1] 5.121677

resultsOutSharp$Rsquared

[1] 0.6189214

# store results
resdHiAnnoyFitAB['Diff sharp', 'RMSE'] <- resultsOutSharp$OOB_RMSE
resdHiAnnoyFitAB['Diff sharp', 'MAE'] <- resultsOutSharp$OOB_MAE
resdHiAnnoyFitAB['Diff sharp', 'Rsquared'] <- resultsOutSharp$Rsquared
resdHiAnnoyPermImpAB$DiffSharp <- resultsOutSharp$conditional_permimp

Plot results

par(mai=c(0,3,0,0))

# plot conditional importance
resultsOutSharp.conimp <- arrange(resultsOutSharp$conditional_permimp, desc(row_number()))

pBar <- ggplot(resultsOutSharp.conimp) + geom_col(aes(x=factor(rownames(resultsOutSharp.conimp), levels=rownames(resultsOutSharp.conimp)), y=CondPermImp), fill=mycolours[2], width=0.5) + labs(x="Variable", y="Conditional variable permutation importance (% highly annoyed)") + ggtitle("dSharpness") + theme(text = element_text(family = "serif"), panel.grid=element_line(color = rgb(235, 235, 235, 100, maxColorValue = 255), linewidth = 0.25, linetype = 2)) + coord_flip()
pBar

if (saveplots){
  ggsave(filename="PtsABdHiAnnoydSharpConPermimp.svg", width=8, height=2.6, path=file.path(outFigPath, "svg"))
  unlink("PtsABdHiAnnoydSharpConPermimp.svg")
  
  ggsave(filename="PtsABdHiAnnoydSharpConPermimp.pdf", width=8, height=2.6, path=file.path(outFigPath, "pdf"))
  unlink("PtsABdHiAnnoydSharpConPermimp.pdf")
}

Selected metric

dSharpVar <- "dSharpAurSHMPowAvgBin"

dTonal loudness and dtonality

Set variables

iVars <- c(diffVar, eventVar, "dTonalECMAAvgMaxLR", "dTonalSHMInt05ExMaxLR", "dTonalSHMIntAvgMaxLR", "dTonalECMA05ExMaxLR", "dTonalAwSHMAvgMaxLR",  "dTonalAwSHM05ExMaxLR", "dTonalAwSHMIntAvgMaxLR",   "dTonalAwSHMInt05ExMaxLR", "dTonLdECMAPowAvgBin", "dTonLdECMA05ExBin", "dTonShpAurSHMPowAvgBin",
           "dTonShpAurSHM05ExBin", "PartTonLdSHMPowAvgBin")
dVar <- "dHighAnnoyPc"

Hyperparameter tuning

p <- mtryTune(dataIn=stimDataNum, iVars=iVars, dVar=dVar, seeds=seeds[1:2],
             ntrees=ntrees, minsplit=minsplit, minbucket=minbucket)

p

seeds <- c(561684, 104798, 1536, 48, 48561)

Selected hyperparameters

ntree <- 1001
mtry <- as.integer(length(iVars)/1.5)

Run model

Train preliminary model

# Tonality with tonal loudness

nperm <- 5

resultsOutTonal1 <- crfReg(dataIn=stimDataNum, iVars=iVars, dVar=dVar, seeds=seeds[1:2], ntree=ntree, mtry=mtry, permImpCondThres=permImpCondThres, nperm=nperm, minsplit=minsplit, minbucket=minbucket)

# print model prediction results
resultsOutTonal1$OOB_RMSE

[1] 7.47396

resultsOutTonal1$OOB_MAE

[1] 5.854586

resultsOutTonal1$Rsquared

[1] 0.5154689

Train multiple seeds model

# Tonality with tonal loudness

resultsOutTonal1 <- multi_crfReg(dataIn=stimDataNum, iVars=iVars, dVar=dVar, seeds=seeds, ntree=ntree, mtry=mtry, permImpCondThres=permImpCondThres, nperm=nperm, minsplit=minsplit, minbucket=minbucket)

# print model prediction results
resultsOutTonal1$OOB_RMSE

[1] 7.473177

resultsOutTonal1$OOB_MAE

[1] 5.855497

resultsOutTonal1$Rsquared

[1] 0.5152894

# store results
resdHiAnnoyFitAB['Diff tonal inc loud', 'RMSE'] <- resultsOutTonal1$OOB_RMSE
resdHiAnnoyFitAB['Diff tonal inc loud', 'MAE'] <- resultsOutTonal1$OOB_MAE
resdHiAnnoyFitAB['Diff tonal inc loud', 'Rsquared'] <- resultsOutTonal1$Rsquared
resdHiAnnoyPermImpAB$DiffTonal1 <- resultsOutTonal1$conditional_permimp

Plot results

par(mai=c(0,3,0,0))

# plot conditional importance
resultsOutTonal1.conimp <- arrange(resultsOutTonal1$conditional_permimp, desc(row_number()))

pBar <- ggplot(resultsOutTonal1.conimp) + geom_col(aes(x=factor(rownames(resultsOutTonal1.conimp), levels=rownames(resultsOutTonal1.conimp)), y=CondPermImp), fill=mycolours[3], width=0.5) + labs(x="Variable", y="Conditional variable permutation importance (% highly annoyed)") + ggtitle("dTonality inc. dtonal loudness") + theme(text = element_text(family = "serif"), panel.grid=element_line(color = rgb(235, 235, 235, 100, maxColorValue = 255), linewidth = 0.25, linetype = 2)) + coord_flip(ylim=c(0, 50))
pBar

if (saveplots){
  ggsave(filename="PtsABdHiAnnoydTonalLdConPermimp.svg", width=8, height=3.6, path=file.path(outFigPath, "svg"))
  unlink("PtsABdHiAnnoydTonalLdConPermimp.svg")
  
  ggsave(filename="PtsABdHiAnnoydTonalLdConPermimp.pdf", width=8, height=3.6, path=file.path(outFigPath, "pdf"))
  unlink("PtsABdHiAnnoydTonalLdConPermimp.pdf")
}

Selected metric

dTonLdVar <- "dTonLdECMAPowAvgBin"

dTonality without dtonal loudness

Set variables

iVars <- c(diffVar, eventVar, "dTonalECMAAvgMaxLR", "dTonalSHMInt05ExMaxLR", "dTonalSHMIntAvgMaxLR", "dTonalECMA05ExMaxLR", "dTonalAwSHMAvgMaxLR",  "dTonalAwSHM05ExMaxLR", "dTonalAwSHMIntAvgMaxLR",   "dTonalAwSHMInt05ExMaxLR")
dVar <- "dHighAnnoyPc"

seeds <- c(410865, 2954, 70812, 203, 7984)

Hyperparameter tuning

p <- mtryTune(dataIn=stimDataNum, iVars=iVars, dVar=dVar, seeds=seeds[1:2],
             ntrees=ntrees, minsplit=minsplit, minbucket=minbucket)

p

NA
NA

Selected hyperparameters

ntree <- 5501
mtry <- as.integer(length(iVars)/1.25)

Run model

Train preliminary model

# Tonality

nperm <- 5

resultsOutTonal2 <- crfReg(dataIn=stimDataNum, iVars=iVars, dVar=dVar, seeds=seeds[1:2], ntree=ntree, mtry=mtry, permImpCondThres=permImpCondThres, nperm=nperm, minsplit=minsplit, minbucket=minbucket)

# print model prediction results
resultsOutTonal2$OOB_RMSE

[1] 7.938544

resultsOutTonal2$OOB_MAE

[1] 6.105294

resultsOutTonal2$Rsquared

[1] 0.4510328

Train multiple seeds model

# Tonality

resultsOutTonal2 <- multi_crfReg(dataIn=stimDataNum, iVars=iVars, dVar=dVar, seeds=seeds, ntree=ntree, mtry=mtry, permImpCondThres=permImpCondThres, nperm=nperm, minsplit=minsplit, minbucket=minbucket)

# print model prediction results
resultsOutTonal2$OOB_RMSE

[1] 7.917518

resultsOutTonal2$OOB_MAE

[1] 6.092962

resultsOutTonal2$Rsquared

[1] 0.4539312

# store results
resdHiAnnoyFitAB['Diff tonal no loud', 'RMSE'] <- resultsOutTonal2$OOB_RMSE
resdHiAnnoyFitAB['Diff tonal no loud', 'MAE'] <- resultsOutTonal2$OOB_MAE
resdHiAnnoyFitAB['Diff tonal no loud', 'Rsquared'] <- resultsOutTonal2$Rsquared
resdHiAnnoyPermImpAB$DiffTonal2 <- resultsOutTonal2$conditional_permimp

Plot results

par(mai=c(0,3,0,0))

# plot conditional importance
resultsOutTonal2.conimp <- arrange(resultsOutTonal2$conditional_permimp, desc(row_number()))

pBar <- ggplot(resultsOutTonal2.conimp) + geom_col(aes(x=factor(rownames(resultsOutTonal2.conimp), levels=rownames(resultsOutTonal2.conimp)), y=CondPermImp), fill=mycolours[3], width=0.5) + labs(x="Variable", y="Conditional variable permutation importance (% highly annoyed)") + ggtitle("dTonality w/o tonal loudness") + theme(text = element_text(family = "serif"), panel.grid=element_line(color = rgb(235, 235, 235, 100, maxColorValue = 255), linewidth = 0.25, linetype = 2)) + coord_flip(ylim=c(0, 50))
pBar

if (saveplots){
  ggsave(filename="PtsABdHiAnnoydTonalConPermimp.svg", width=8, height=2.8, path=file.path(outFigPath, "svg"))
  unlink("PtsABdHiAnnoydTonalConPermimp.svg")
  
  ggsave(filename="PtsABdHiAnnoydTonalConPermimp.pdf", width=8, height=2.8, path=file.path(outFigPath, "pdf"))
  unlink("PtsABdHiAnnoydTonalConPermimp.pdf")
}

Selected metric

dTonalVar <- "dTonalAwSHMIntAvgMaxLR"

dFluctuation strength

Set variables

# Fluctuation strength
iVars <- c(diffVar, eventVar, "dFluctECMA10ExBin", "dFluctECMA05ExBin", "dFluctOV10ExMaxLR", "dFluctOV05ExMaxLR")
dVar <- "dHighAnnoyPc"

seeds <- c(418657, 84, 1630, 18659, 3687)

Hyperparameter tuning

p <- mtryTune(dataIn=stimDataNum, iVars=iVars, dVar=dVar, seeds=seeds[1:2],
             ntrees=ntrees, minsplit=minsplit, minbucket=minbucket)

p

Selected hyperparameters

ntree <- 2501
mtry <- 3

Run model

Train preliminary model

nperm <- 5

resultsOutFluct <- crfReg(dataIn=stimDataNum, iVars=iVars, dVar=dVar, seeds=seeds[1:2], ntree=ntree, mtry=mtry, permImpCondThres=permImpCondThres, nperm=nperm, minsplit=minsplit, minbucket=minbucket)

# print model prediction results
resultsOutFluct$OOB_RMSE

[1] 7.887561

resultsOutFluct$OOB_MAE

[1] 6.188833

resultsOutFluct$Rsquared

[1] 0.4586595

Train multiple seeds model

resultsOutFluct <- multi_crfReg(dataIn=stimDataNum, iVars=iVars, dVar=dVar, seeds=seeds, ntree=ntree, mtry=mtry, permImpCondThres=permImpCondThres, nperm=nperm, minsplit=minsplit, minbucket=minbucket)

# print model prediction results
resultsOutFluct$OOB_RMSE

[1] 7.875212

resultsOutFluct$OOB_MAE

[1] 6.181713

resultsOutFluct$Rsquared

[1] 0.4605064

# store results
resdHiAnnoyFitAB['Diff fluct', 'RMSE'] <- resultsOutFluct$OOB_RMSE
resdHiAnnoyFitAB['Diff fluct', 'MAE'] <- resultsOutFluct$OOB_MAE
resdHiAnnoyFitAB['Diff fluct', 'Rsquared'] <- resultsOutFluct$Rsquared
resdHiAnnoyPermImpAB$DiffFluct <- resultsOutFluct$conditional_permimp

Plot results

par(mai=c(0,3,0,0))

# plot conditional importance
resultsOutFluct.conimp <- arrange(resultsOutFluct$conditional_permimp, desc(row_number()))

pBar <- ggplot(resultsOutFluct.conimp) + geom_col(aes(x=factor(rownames(resultsOutFluct.conimp), levels=rownames(resultsOutFluct.conimp)), y=CondPermImp), fill=mycolours[4], width=0.5) + labs(x="Variable", y="Conditional variable permutation importance (% highly annoyed)") + ggtitle("dFluctuation strength") + theme(text = element_text(family = "serif"), panel.grid=element_line(color = rgb(235, 235, 235, 100, maxColorValue = 255), linewidth = 0.25, linetype = 2)) + coord_flip()
pBar

if (saveplots){
  ggsave(filename="PtsABdHiAnnoydFluctConPermimp.svg", width=8, height=2, path=file.path(outFigPath, "svg"))
  unlink("PtsABdHiAnnoyFluctConPermimp.svg")
  
  ggsave(filename="PtsABdHiAnnoydFluctConPermimp.pdf", width=8, height=2, path=file.path(outFigPath, "pdf"))
  unlink("PtsABdHiAnnoyFluctConPermimp.pdf")
}

Selected metric

dFluctVar <- "dFluctOV10ExMaxLR"

dRoughness

Set variables

# Roughness
iVars <- c(diffVar, eventVar, "dRoughECMA10ExBin", "dRoughECMA05ExBin", "dRoughFZ10ExMaxLR", "dRoughFZ05ExMaxLR")
dVar <- "dHighAnnoyPc"

seeds <- c(69851, 85109, 410986, 1563, 896)

Hyperparameter tuning

p <- mtryTune(dataIn=stimDataNum, iVars=iVars, dVar=dVar, seeds=seeds[1:2],
             ntrees=ntrees, minsplit=minsplit, minbucket=minbucket)

p

Selected hyperparameters

ntree <- 501
mtry <- 3

Run model

Train preliminary model

nperm <- 5

resultsOutRough <- crfReg(dataIn=stimDataNum, iVars=iVars, dVar=dVar, seeds=seeds[1:2], ntree=ntree, mtry=mtry, permImpCondThres=permImpCondThres, nperm=nperm, minsplit=minsplit, minbucket=minbucket)

# print model prediction results
resultsOutRough$OOB_RMSE

[1] 7.82689

resultsOutRough$OOB_MAE

[1] 6.109651

resultsOutRough$Rsquared

[1] 0.4701734

Train multiple seeds model

resultsOutRough <- multi_crfReg(dataIn=stimDataNum, iVars=iVars, dVar=dVar, seeds=seeds, ntree=ntree, mtry=mtry, permImpCondThres=permImpCondThres, nperm=nperm, minsplit=minsplit, minbucket=minbucket)

# print model prediction results
resultsOutRough$OOB_RMSE

[1] 7.790608

resultsOutRough$OOB_MAE

[1] 6.100323

resultsOutRough$Rsquared

[1] 0.4760393

# store results
resdHiAnnoyFitAB['Diff rough', 'RMSE'] <- resultsOutRough$OOB_RMSE
resdHiAnnoyFitAB['Diff rough', 'MAE'] <- resultsOutRough$OOB_MAE
resdHiAnnoyFitAB['Diff rough', 'Rsquared'] <- resultsOutRough$Rsquared
resdHiAnnoyPermImpAB$DiffRough <- resultsOutRough$conditional_permimp

Plot results

par(mai=c(0,3,0,0))

# plot conditional importance
resultsOutRough.conimp <- arrange(resultsOutRough$conditional_permimp, desc(row_number()))

pBar <- ggplot(resultsOutRough.conimp) + geom_col(aes(x=factor(rownames(resultsOutRough.conimp), levels=rownames(resultsOutRough.conimp)), y=CondPermImp), fill=mycolours[5], width=0.5) + labs(x="Variable", y="Conditional variable permutation importance (% highly annoyed)") + ggtitle("dRoughness") + theme(text = element_text(family = "serif"), panel.grid=element_line(color = rgb(235, 235, 235, 100, maxColorValue = 255), linewidth = 0.25, linetype = 2)) + coord_flip()
pBar

if (saveplots){
  ggsave(filename="PtsABdHiAnnoydRoughConPermimp.svg", width=8, height=2, path=file.path(outFigPath, "svg"))
  unlink("PtsABdHiAnnoydRoughConPermimp.svg")
  
  ggsave(filename="PtsABdHiAnnoydRoughConPermimp.pdf", width=8, height=2, path=file.path(outFigPath, "pdf"))
  unlink("PtsABdHiAnnoydRoughConPermimp.pdf")
}

Selected metric

dRoughVar <- "dRoughFZ05ExMaxLR"

dImpulsiveness

Set variables

# Impulsiveness
iVars <- c(diffVar, eventVar, "dImpulsSHMAvgMaxLR", "dImpulsSHM05ExMaxLR", "dImpulsSHMPowAvgMaxLR",
           "dImpulsLoudWZAvgMaxLR", "dImpulsLoudWZ05ExMaxLR", "dImpulsLoudWZPowAvgMaxLR",
           "dImpulsLoudWECMAAvgBin", "dImpulsLoudWECMA05ExBin", "dImpulsLoudWECMAPowAvgBin")
dVar <- "dHighAnnoyPc"

seeds <- c(418659, 7805, 38475, 65834, 1653)

Hyperparameter tuning

p <- mtryTune(dataIn=stimDataNum, iVars=iVars, dVar=dVar, seeds=seeds[1:2],
             ntrees=ntrees, minsplit=minsplit, minbucket=minbucket)

p

Selected hyperparameters

ntree <- 5501
mtry <- as.integer(length(iVars)/1.25)

Run model

Train preliminary model

nperm <- 5

resultsOutImpuls <- crfReg(dataIn=stimDataNum, iVars=iVars, dVar=dVar, seeds=seeds[1:2], ntree=ntree, mtry=mtry, permImpCondThres=permImpCondThres, nperm=nperm, minsplit=minsplit, minbucket=minbucket)

# print model prediction results
resultsOutImpuls$OOB_RMSE

[1] 7.444946

resultsOutImpuls$OOB_MAE

[1] 5.833038

resultsOutImpuls$Rsquared

[1] 0.5179816

Train multiple seeds model

resultsOutImpuls <- multi_crfReg(dataIn=stimDataNum, iVars=iVars, dVar=dVar, seeds=seeds, ntree=ntree, mtry=mtry, permImpCondThres=permImpCondThres, nperm=nperm, minsplit=minsplit, minbucket=minbucket)

# print model prediction results
resultsOutImpuls$OOB_RMSE

[1] 7.439293

resultsOutImpuls$OOB_MAE

[1] 5.831675

resultsOutImpuls$Rsquared

[1] 0.5186948

# store results
resdHiAnnoyFitAB['Diff impuls', 'RMSE'] <- resultsOutImpuls$OOB_RMSE
resdHiAnnoyFitAB['Diff impuls', 'MAE'] <- resultsOutImpuls$OOB_MAE
resdHiAnnoyFitAB['Diff impuls', 'Rsquared'] <- resultsOutImpuls$Rsquared
resdHiAnnoyPermImpAB$DiffImpuls <- resultsOutImpuls$conditional_permimp

Plot results

par(mai=c(0,3,0,0))

# plot conditional importance
resultsOutImpuls.conimp <- arrange(resultsOutImpuls$conditional_permimp, desc(row_number()))

pBar <- ggplot(resultsOutImpuls.conimp) + geom_col(aes(x=factor(rownames(resultsOutImpuls.conimp), levels=rownames(resultsOutImpuls.conimp)), y=CondPermImp), fill=mycolours[6], width=0.5) + labs(x="Variable", y="Conditional variable permutation importance (% highly annoyed)") + ggtitle("dImpulsiveness") + theme(text = element_text(family = "serif"), panel.grid=element_line(color = rgb(235, 235, 235, 100, maxColorValue = 255), linewidth = 0.25, linetype = 2)) + coord_flip()
pBar

if (saveplots){
  ggsave(filename="PtsABdHiAnnoydImpulsConPermimp.svg", width=8, height=3.6, path=file.path(outFigPath, "svg"))
  unlink("PtsABdHiAnnoydImpulsConPermimp.svg")
  
  ggsave(filename="PtsABdHiAnnoydImpulsConPermimp.pdf", width=8, height=3.6, path=file.path(outFigPath, "pdf"))
  unlink("PtsABdHiAnnoydImpulsConPermimp.pdf")
}

Selected metric

dImpulsVar <- "dImpulsLoudWZAvgMaxLR"

dSQM and loudness comparison

Now the highest importance dSQMs are ranked against each other, controlling for loudness difference.

Include dtonal loudness

Set variables

iVars <- c(diffVar, eventVar, dSharpVar, dTonLdVar, dFluctVar, dRoughVar, dImpulsVar)
dVar <- "dHighAnnoyPc"

seeds <- c(98465, 54163, 6541, 36485, 849675)

Hyperparameter tuning

p <- mtryTune(dataIn=stimDataNum, iVars=iVars, dVar=dVar, seeds=seeds[1:2],
             ntrees=ntrees, minsplit=minsplit, minbucket=minbucket)

p

Selected hyperparameters

ntree <- 251
mtry <- as.integer(length(iVars)/1.75)

Run model

Train preliminary model

nperm <- 5

resultsOutSQMs1 <- crfReg(dataIn=stimDataNum, iVars=iVars, dVar=dVar, seeds=seeds[1:2], ntree=ntree, mtry=mtry, permImpCondThres=permImpCondThres, nperm=nperm, minsplit=minsplit, minbucket=minbucket)

# print model prediction results
resultsOutSQMs1$OOB_RMSE

[1] 7.123136

resultsOutSQMs1$OOB_MAE

[1] 5.628937

resultsOutSQMs1$Rsquared

[1] 0.5583272

Train multiple seeds model

resultsOutSQMs1 <- multi_crfReg(dataIn=stimDataNum, iVars=iVars, dVar=dVar, seeds=seeds, ntree=ntree, mtry=mtry, permImpCondThres=permImpCondThres, nperm=nperm, minsplit=minsplit, minbucket=minbucket)

# print model prediction results
resultsOutSQMs1$OOB_RMSE

[1] 7.100365

resultsOutSQMs1$OOB_MAE

[1] 5.589069

resultsOutSQMs1$Rsquared

[1] 0.5618113

# store results
resdHiAnnoyFitAB['Diff SQMs inc tonal loud', 'RMSE'] <- resultsOutSQMs1$OOB_RMSE
resdHiAnnoyFitAB['Diff SQMs inc tonal loud', 'MAE'] <- resultsOutSQMs1$OOB_MAE
resdHiAnnoyFitAB['Diff SQMs inc tonal loud', 'Rsquared'] <- resultsOutSQMs1$Rsquared
resdHiAnnoyPermImpAB$DiffSQMs1 <- resultsOutSQMs1$conditional_permimp

Plot results

par(mai=c(0,3,0,0))

# plot conditional importance
resultsOutSQMs1.conimp <- arrange(resultsOutSQMs1$conditional_permimp, desc(row_number()))

pBar <- ggplot(resultsOutSQMs1.conimp) + geom_col(aes(x=factor(rownames(resultsOutSQMs1.conimp), levels=rownames(resultsOutSQMs1.conimp)), y=CondPermImp), fill=mycolours[7], width=0.5) + labs(x="Variable", y="Conditional variable permutation importance (% highly annoyed)") + theme(text = element_text(family = "serif"), panel.grid=element_line(color = rgb(235, 235, 235, 100, maxColorValue = 255), linewidth = 0.25, linetype = 2)) + coord_flip(ylim=c(0, 50))
pBar

if (saveplots){
  ggsave(filename="PtsABdHiAnnoyDiffSQMsTonLdConPermimp.svg", width=8, height=2.4, path=file.path(outFigPath, "svg"))
  unlink("PtsABdHiAnnoyDiffSQMsTonLdConPermimp.svg")
  
  ggsave(filename="PtsABdHiAnnoyDiffSQMsTonLdConPermimp.pdf", width=8, height=2.4, path=file.path(outFigPath, "pdf"))
  unlink("PtsABdHiAnnoyDiffSQMsTonLdConPermimp.pdf")
}

Exclude tonal loudness

Set variables

iVars <- c(diffVar, eventVar, dSharpVar, dTonalVar, dFluctVar, dRoughVar, dImpulsVar)
dVar <- "dHighAnnoyPc"

seeds <- c(49865, 7852, 845961, 410583, 36748)

Hyperparameter tuning

p <- mtryTune(dataIn=stimDataNum, iVars=iVars, dVar=dVar, seeds=seeds[1:2],
             ntrees=ntrees, minsplit=minsplit, minbucket=minbucket)

p

Selected hyperparameters

ntree <- 501
mtry <- as.integer(length(iVars)/1.75)

Run model

Train preliminary model

nperm <- 5

resultsOutSQMs2 <- crfReg(dataIn=stimDataNum, iVars=iVars, dVar=dVar, seeds=seeds[1:2], ntree=ntree, mtry=mtry, permImpCondThres=permImpCondThres, nperm=nperm, minsplit=minsplit, minbucket=minbucket)

# print model prediction results
resultsOutSQMs2$OOB_RMSE

[1] 7.09875

resultsOutSQMs2$OOB_MAE

[1] 5.575011

resultsOutSQMs2$Rsquared

[1] 0.5623681

Train multiple seeds model

resultsOutSQMs2 <- multi_crfReg(dataIn=stimDataNum, iVars=iVars, dVar=dVar, seeds=seeds, ntree=ntree, mtry=mtry, permImpCondThres=permImpCondThres, nperm=nperm, minsplit=minsplit, minbucket=minbucket)

# print model prediction results
resultsOutSQMs2$OOB_RMSE

[1] 7.070765

resultsOutSQMs2$OOB_MAE

[1] 5.560362

resultsOutSQMs2$Rsquared

[1] 0.56598

# store results
resdHiAnnoyFitAB['Diff SQMs no tonal loud', 'RMSE'] <- resultsOutSQMs2$OOB_RMSE
resdHiAnnoyFitAB['Diff SQMs no tonal loud', 'MAE'] <- resultsOutSQMs2$OOB_MAE
resdHiAnnoyFitAB['Diff SQMs no tonal loud', 'Rsquared'] <- resultsOutSQMs2$Rsquared
resdHiAnnoyPermImpAB$DiffSQMs2 <- resultsOutSQMs2$conditional_permimp

Plot results

par(mai=c(0,3,0,0))

# plot conditional importance
resultsOutSQMs2.conimp <- arrange(resultsOutSQMs2$conditional_permimp, desc(row_number()))

pBar <- ggplot(resultsOutSQMs2.conimp) + geom_col(aes(x=factor(rownames(resultsOutSQMs2.conimp), levels=rownames(resultsOutSQMs2.conimp)), y=CondPermImp), fill=mycolours[7], width=0.5) + labs(x="Variable", y="Conditional variable permutation importance (% highly annoyed)") + theme(text = element_text(family = "serif"), panel.grid=element_line(color = rgb(235, 235, 235, 100, maxColorValue = 255), linewidth = 0.25, linetype = 2)) + coord_flip(ylim=c(0, 50))
pBar

if (saveplots){
  ggsave(filename="PtsABdHiAnnoyDiffSQMsNoTonLdConPermimp.svg", width=8, height=2.4, path=file.path(outFigPath, "svg"))
  unlink("PtsABdHiAnnoyDiffSQMsNoTonLdConPermimp.svg")
  
  ggsave(filename="PtsABdHiAnnoyDiffSQMsNoTonLdConPermimp.pdf", width=8, height=2.4, path=file.path(outFigPath, "pdf"))
  unlink("PtsABdHiAnnoyDiffSQMsNoTonLdConPermimp.pdf")
}

dPsychoacoustic annoyance metrics

Set variables

iVars <- c(ambVar, "dPsychAnnoyWidmann", "dPsychAnnoyMore", "dPsychAnnoyDi", "dPsychAnnoyTorija", "dPsychAnnoyWillemsen", "dPsychAnnoyBoucher")
dVar <- "dHighAnnoyPc"

seeds <- c(835702, 54, 470912, 652, 55297)

Hyperparameter tuning

p <- mtryTune(dataIn=stimDataNum, iVars=iVars, dVar=dVar, seeds=seeds[1:2],
              ntrees=ntrees, minsplit=minsplit, minbucket=minbucket)

p

Selected hyperparameters

ntree <- 251
mtry <-  4

Run model

Train preliminary model

nperm <- 5

resultsOutPA <- crfReg(dataIn=stimDataNum, iVars=iVars, dVar=dVar, seeds=seeds[1:2], ntree=ntree, mtry=mtry, permImpCondThres=permImpCondThres, nperm=nperm, minsplit=minsplit, minbucket=minbucket)

# print model prediction results
resultsOutPA$OOB_RMSE

[1] 7.901164

resultsOutPA$OOB_MAE

[1] 5.982122

resultsOutPA$Rsquared

[1] 0.4558801

Train multiple seeds model

resultsOutPA <- multi_crfReg(dataIn=stimDataNum, iVars=iVars, dVar=dVar, seeds=seeds, ntree=ntree, mtry=mtry, permImpCondThres=permImpCondThres, nperm=nperm, minsplit=minsplit, minbucket=minbucket)

# print model prediction results
resultsOutPA$OOB_RMSE

[1] 7.89992

resultsOutPA$OOB_MAE

[1] 5.953038

resultsOutPA$Rsquared

[1] 0.4563591

# store results
resdHiAnnoyFitAB['Diff Psychoacoustic annoyance', 'RMSE'] <- resultsOutPA$OOB_RMSE
resdHiAnnoyFitAB['Diff Psychoacoustic annoyance', 'MAE'] <- resultsOutPA$OOB_MAE
resdHiAnnoyFitAB['Diff Psychoacoustic annoyance', 'Rsquared'] <- resultsOutPA$Rsquared
resdHiAnnoyPermImpAB$AbsPA <- resultsOutPA$conditional_permimp

Plot results

par(mai=c(0,3,0,0))

# plot conditional importance
resultsOutPA.conimp <- arrange(resultsOutPA$conditional_permimp, desc(row_number()))

pBar <- ggplot(resultsOutPA.conimp) + geom_col(aes(x=factor(rownames(resultsOutPA.conimp), levels=rownames(resultsOutPA.conimp)), y=CondPermImp), fill=mycolours[10], width=0.5) + labs(x="Variable", y="Conditional variable permutation importance (% highly annoyed)") + theme(text = element_text(family = "serif"), panel.grid=element_line(color = rgb(235, 235, 235, 100, maxColorValue = 255), linewidth = 0.25, linetype = 2)) + coord_flip(ylim=c(0, 50))
pBar

if (saveplots){
  ggsave(filename="PtsABdHiAnnoydPAConPermimp.svg", width=8, height=2.4, path=file.path(outFigPath, "svg"))
  unlink("PtsABdHiAnnoydPAConPermimp.svg")
  
  ggsave(filename="PtsABdHiAnnoydPAConPermimp.pdf", width=8, height=2.4, path=file.path(outFigPath, "pdf"))
  unlink("PtsABdHiAnnoydPAConPermimp.pdf")
}

Save the results outputs to file

if (savedata){
  utils::write.csv(resdHiAnnoyFitAB, paste(outDataPath, "\\PtsABCRFdHiAnnoyOOBFit.csv", sep=""))
  ii <- 0
  temp = list()
  for (res in resdHiAnnoyPermImpAB){
    ii <- ii + 1
    temp[[ii]] <- as.data.frame(resdHiAnnoyPermImpAB[ii])
    names(temp[[ii]]) <- names(resdHiAnnoyPermImpAB[ii])
  }
  openxlsx::write.xlsx(temp, paste(outDataPath, "\\PtsABCRFdHiAnnoyConPermimp.xlsx",
                                   sep=""),
                       rowNames=TRUE)
}

Parts A&B summary

Summary of results for Parts A & B combined

With tonal loudness

Absolute variables

# combine the annoyance perm importance results

# convert each result to a tibble with rownames added to a column, renaming the data column to 'dAnnoy' etc.
resdAnnoyMnAbsPermImpTblAB <- as.data.frame(resdAnnoyMnPermImpAB$AbsSQMs1/max(resdAnnoyMnPermImpAB$AbsSQMs1)) |>
  tibble::rownames_to_column(var='Variable')
colnames(resdAnnoyMnAbsPermImpTblAB)[2] <- "dAnnoy"

resdHiAnnoyAbsPermImpTblAB <- as.data.frame(resdHiAnnoyPermImpAB$AbsSQMs1/max(resdHiAnnoyPermImpAB$AbsSQMs1)) |>
  tibble::rownames_to_column(var='Variable')
colnames(resdHiAnnoyAbsPermImpTblAB)[2] <- "dHiAnnoy"

# merge the dataframes
resAbsPermImpTblAB <- list(resdAnnoyMnAbsPermImpTblAB, resdHiAnnoyAbsPermImpTblAB) |>
  purrr::reduce(merge, by = c('Variable'), all = T)

# rename the columns
colnames(resAbsPermImpTblAB)[2:3] <- c("Mean change in annoyance", "%HA | HA' (amb.)")
resAbsPermImpTblAB[is.na(resAbsPermImpTblAB)] <- 0

resAbsAB <- tidyr::pivot_longer(resAbsPermImpTblAB, cols=-Variable, names_to="Outcome", values_to="Imp")

# reorder res tibble, descending by the variable Imp grouped sum and create column with new group order as a factor
resAbsAB <- resAbsAB |> mutate(Variable_sum = sum(Imp), .by=Variable) |> arrange(desc(Variable_sum)) |> group_by(Variable_sum, Variable) |>
   mutate(Order = cur_group_id()) |> mutate(Order = as.factor(Order)) |> arrange(desc(Order))

# Reorder outcome levels
resAbsAB$Outcome <- factor(resAbsAB$Outcome, levels=c("Mean change in annoyance", "%HA | HA' (amb.)"))

# plot res as horizontal bar chart, with Imp as y axis, Variable as x axis, Outcome as fill, and Variable_sum as order, relabel x axis with Variable names
pBar <- ggplot(resAbsAB) + geom_col(aes(fill=Outcome, y=Imp, x=Order), colour='grey35', linewidth=0,  width=0.75, show.legend=TRUE) + labs(x="Variable", y="Normalised conditional variable\npermutation importance") + theme(text = element_text(family = "serif"), panel.grid=element_line(color = rgb(235, 235, 235, 100, maxColorValue = 255), linewidth = 0.25, linetype = 2)) + coord_flip(ylim=c(0, 1)) + scale_fill_manual(values=mycolours, labels=c(expression(paste(bar(Delta~A))), "%HA | HA' (amb.)")) + scale_x_discrete(labels=unique(rev(resAbsAB$Variable))) + guides(fill=guide_legend(title='Outcome'))
pBar + scale_y_continuous(breaks=seq(0, 1, by=0.5))

if (saveplots){
  ggsave(filename="PtsABcrfAbsSQMsSummary.svg", width=8, height=4, path=file.path(outFigPath, "svg"))
  unlink("PtsABcrfAbsSQMsSummary.svg")

  ggsave(filename="PtsABcrfAbsSQMsSummary.pdf", width=8, height=4, path=file.path(outFigPath, "pdf"))
  unlink("PtsABcrfAbsSQMsSummary.pdf")
}


# plot res as horizontal bar chart, with Imp as y axis, Variable as x axis, Outcome as fill, and Variable_sum as order, relabel x axis with Variable names
pBar <- ggplot(resAbsAB) + geom_col(aes(fill=Outcome, y=Imp, x=Order), colour='grey35', linewidth=0,  width=0.75, show.legend=TRUE) + labs(x="Variable", y="Normalised conditional variable\npermutation importance") + theme(text = element_text(family = "serif"), panel.grid=element_line(color = rgb(235, 235, 235, 100, maxColorValue = 255), linewidth = 0.25, linetype = 2), legend.position = "top") + coord_flip(ylim=c(0, 1)) + scale_fill_manual(values=mycolours, labels=c(expression(paste(bar(Delta~A))), "%HA | HA' (amb.)")) + scale_x_discrete(labels=unique(rev(resAbsAB$Variable))) + guides(fill=guide_legend(title='Outcome', nrow=2, ncol=1))
pBar + scale_y_continuous(breaks=seq(0, 1, by=0.5))

if (saveplots){
  ggsave(filename="PtsABcrfAbsSQMsSummaryNw.svg", width=4, height=4, path=file.path(outFigPath, "svg"))
  unlink("PtsABcrfAbsSQMsSummary.svg")

  ggsave(filename="PtsABcrfAbsSQMsSummaryNw.pdf", width=4, height=4, path=file.path(outFigPath, "pdf"))
  unlink("PtsABcrfAbsSQMsSummary.pdf")
}

Difference variables

# combine the annoyance perm importance results

# convert each result to a tibble with rownames added to a column, renaming the data column to 'dAnnoy' etc.
resdAnnoyMnDiffPermImpTblAB <- as.data.frame(resdAnnoyMnPermImpAB$DiffSQMs1/max(resdAnnoyMnPermImpAB$DiffSQMs1)) |>
  tibble::rownames_to_column(var='Variable')
colnames(resdAnnoyMnDiffPermImpTblAB)[2] <- "dAnnoy"

resdHiAnnoyDiffPermImpTblAB <- as.data.frame(resdHiAnnoyPermImpAB$DiffSQMs1/max(resdHiAnnoyPermImpAB$DiffSQMs1)) |>
  tibble::rownames_to_column(var='Variable')
colnames(resdHiAnnoyDiffPermImpTblAB)[2] <- "dHiAnnoy"

# merge the dataframes
resDiffPermImpTblAB <- list(resdAnnoyMnDiffPermImpTblAB, resdHiAnnoyDiffPermImpTblAB) |>
  purrr::reduce(merge, by = c('Variable'), all = T)

# rename the columns
colnames(resDiffPermImpTblAB)[2:3] <- c("Mean change in annoyance", "%HA | HA' (amb.)")
resDiffPermImpTblAB[is.na(resDiffPermImpTblAB)] <- 0

resDiffAB <- tidyr::pivot_longer(resDiffPermImpTblAB, cols=-Variable, names_to="Outcome", values_to="Imp")

# reorder res tibble, descending by the variable Imp grouped sum and create column with new group order as a factor
resDiffAB <- resDiffAB |> mutate(Variable_sum = sum(Imp), .by=Variable) |> arrange(desc(Variable_sum)) |> group_by(Variable_sum, Variable) |>
   mutate(Order = cur_group_id()) |> mutate(Order = as.factor(Order)) |> arrange(desc(Order))

# Reorder outcome levels
resDiffAB$Outcome <- factor(resDiffAB$Outcome, levels=c("Mean change in annoyance", "%HA | HA' (amb.)"))

# plot res as horizontal bar chart, with Imp as y axis, Variable as x axis, Outcome as fill, and Variable_sum as order, relabel x axis with Variable names
pBar <- ggplot(resDiffAB) + geom_col(aes(fill=Outcome, y=Imp, x=Order), colour='grey35', linewidth=0,  width=0.75, show.legend=TRUE) + labs(x="Variable", y="Normalised conditional variable permutation importance") + theme(text = element_text(family = "serif"), panel.grid=element_line(color = rgb(235, 235, 235, 100, maxColorValue = 255), linewidth = 0.25, linetype = 2), legend.position = "right") + coord_flip(ylim=c(-0.1, 1.1)) + scale_fill_manual(values=mycolours) + scale_x_discrete(labels=unique(rev(resDiffAB$Variable)))
pBar + scale_y_continuous(breaks=seq(0, 1, by=0.5))

if (saveplots){
  ggsave(filename="PtsABcrfDiffSQMsSummary.svg", width=8, height=3, path=file.path(outFigPath, "svg"))
  unlink("PtsABcrfDiffSQMsSummary.svg")
  
  ggsave(filename="PtsABcrfDiffSQMsSummary.pdf", width=8, height=3, path=file.path(outFigPath, "pdf"))
  unlink("PtsABcrfDiffSQMsSummary.pdf")
}

No tonal loudness

Absolute variables

# combine the annoyance perm importance results

# convert each result to a tibble with rownames added to a column, renaming the data column to 'dAnnoy' etc.
resdAnnoyMnAbsPermImpNoTonLdTblAB <- as.data.frame(resdAnnoyMnPermImpAB$AbsSQMs2/max(resdAnnoyMnPermImpAB$AbsSQMs2)) |>
  tibble::rownames_to_column(var='Variable')
colnames(resdAnnoyMnAbsPermImpNoTonLdTblAB)[2] <- "dAnnoy"

resdHiAnnoyAbsPermImpNoTonLdTblAB <- as.data.frame(resdHiAnnoyPermImpAB$AbsSQMs2/max(resdHiAnnoyPermImpAB$AbsSQMs2)) |>
  tibble::rownames_to_column(var='Variable')
colnames(resdHiAnnoyAbsPermImpNoTonLdTblAB)[2] <- "dHiAnnoy"

# merge the dataframes
resAbsPermImpNoTonLdTblAB <- list(resdAnnoyMnAbsPermImpNoTonLdTblAB, resdHiAnnoyAbsPermImpNoTonLdTblAB) |>
  purrr::reduce(merge, by = c('Variable'), all = T)

# rename the columns
colnames(resAbsPermImpNoTonLdTblAB)[2:3] <- c("Mean change in annoyance", "%HA | HA' (amb.)")
resAbsPermImpNoTonLdTblAB[is.na(resAbsPermImpNoTonLdTblAB)] <- 0

resAbsNoTonLdAB <- tidyr::pivot_longer(resAbsPermImpNoTonLdTblAB, cols=-Variable, names_to="Outcome", values_to="Imp")

# reorder res tibble, descending by the variable Imp grouped sum and create column with new group order as a factor
resAbsNoTonLdAB <- resAbsNoTonLdAB |> mutate(Variable_sum = sum(Imp), .by=Variable) |> arrange(desc(Variable_sum)) |> group_by(Variable_sum, Variable) |>
   mutate(Order = cur_group_id()) |> mutate(Order = as.factor(Order)) |> arrange(desc(Order))

# Reorder outcome levels
resAbsNoTonLdAB$Outcome <- factor(resAbsNoTonLdAB$Outcome, levels=c("Mean change in annoyance", "%HA | HA' (amb.)"))

# plot res as horizontal bar chart, with Imp as y axis, Variable as x axis, Outcome as fill, and Variable_sum as order, relabel x axis with Variable names
pBar <- ggplot(resAbsNoTonLdAB) + geom_col(aes(fill=Outcome, y=Imp, x=Order), colour='grey35', linewidth=0,  width=0.75, show.legend=TRUE) + labs(x="Variable", y="Normalised conditional variable permutation importance") + theme(text = element_text(family = "serif"), panel.grid=element_line(color = rgb(235, 235, 235, 100, maxColorValue = 255), linewidth = 0.25, linetype = 2), legend.position = "right") + coord_flip(ylim=c(-0.1, 1.1)) + scale_fill_manual(values=mycolours) + scale_x_discrete(labels=unique(rev(resAbsNoTonLdAB$Variable)))
pBar + scale_y_continuous(breaks=seq(0, 1, by=0.5))

if (saveplots){
  ggsave(filename="PtsABcrfAbsSQMsNoTonLdSummary.svg", width=8, height=3, path=file.path(outFigPath, "svg"))
  unlink("PtsABcrfAbsSQMsNoTonLdSummary.svg")
  
  ggsave(filename="PtsABcrfAbsSQMsNoTonLdSummary.pdf", width=8, height=3, path=file.path(outFigPath, "pdf"))
  unlink("PtsABcrfAbsSQMsNoTonLdSummary.pdf")
}

Difference variables

# combine the annoyance perm importance results

# convert each result to a tibble with rownames added to a column, renaming the data column to 'dAnnoy' etc.
resdAnnoyMnDiffPermImpNoTonLdTblAB <- as.data.frame(resdAnnoyMnPermImpAB$DiffSQMs2/max(resdAnnoyMnPermImpAB$DiffSQMs2)) |>
  tibble::rownames_to_column(var='Variable')
colnames(resdAnnoyMnDiffPermImpNoTonLdTblAB)[2] <- "dAnnoy"

resdHiAnnoyDiffPermImpNoTonLdTblAB <- as.data.frame(resdHiAnnoyPermImpAB$DiffSQMs2/max(resdHiAnnoyPermImpAB$DiffSQMs2)) |>
  tibble::rownames_to_column(var='Variable')
colnames(resdHiAnnoyDiffPermImpNoTonLdTblAB)[2] <- "dHiAnnoy"

# merge the dataframes
resDiffPermImpNoTonLdTblAB <- list(resdAnnoyMnDiffPermImpNoTonLdTblAB, resdHiAnnoyDiffPermImpNoTonLdTblAB) |>
  purrr::reduce(merge, by = c('Variable'), all = T)

# rename the columns
colnames(resDiffPermImpNoTonLdTblAB)[2:3] <- c("Mean change in annoyance", "%HA | HA' (amb.)")
resDiffPermImpNoTonLdTblAB[is.na(resDiffPermImpNoTonLdTblAB)] <- 0

resDiffNoTonLdAB <- tidyr::pivot_longer(resDiffPermImpNoTonLdTblAB, cols=-Variable, names_to="Outcome", values_to="Imp")

# reorder res tibble, descending by the variable Imp grouped sum and create column with new group order as a factor
resDiffNoTonLdAB <- resDiffNoTonLdAB |> mutate(Variable_sum = sum(Imp), .by=Variable) |> arrange(desc(Variable_sum)) |> group_by(Variable_sum, Variable) |>
   mutate(Order = cur_group_id()) |> mutate(Order = as.factor(Order)) |> arrange(desc(Order))

# Reorder outcome levels
resDiffNoTonLdAB$Outcome <- factor(resDiffNoTonLdAB$Outcome, levels=c("Mean change in annoyance", "%HA | HA' (amb.)"))

# plot res as horizontal bar chart, with Imp as y axis, Variable as x axis, Outcome as fill, and Variable_sum as order, relabel x axis with Variable names
pBar <- ggplot(resDiffNoTonLdAB) + geom_col(aes(fill=Outcome, y=Imp, x=Order), colour='grey35', linewidth=0,  width=0.75, show.legend=TRUE) + labs(x="Variable", y="Normalised conditional variable permutation importance") + theme(text = element_text(family = "serif"), panel.grid=element_line(color = rgb(235, 235, 235, 100, maxColorValue = 255), linewidth = 0.25, linetype = 2), legend.position = "right") + coord_flip(ylim=c(-0.1, 1.1)) + scale_fill_manual(values=mycolours) + scale_x_discrete(labels=unique(rev(resDiffNoTonLdAB$Variable)))
pBar + scale_y_continuous(breaks=seq(0, 1, by=0.5))

if (saveplots){
  ggsave(filename="PtsABcrfDiffSQMsNoTonLdSummary.svg", width=8, height=3, path=file.path(outFigPath, "svg"))
  unlink("PtsABcrfDiffSQMsNoTonLdSummary.svg")
  
  ggsave(filename="PtsABcrfDiffSQMsNoTonLdSummary.pdf", width=8, height=3, path=file.path(outFigPath, "pdf"))
  unlink("PtsABcrfDiffSQMsNoTonLdSummary.pdf")
}

Save the results outputs to file

# Make a list of the summary results
resSummary <- list(resAbsAB, resDiffAB, resAbsNoTonLdAB, resDiffNoTonLdAB)

# Save the results
if (savedata){
  ii <- 0
  temp = list()
  for (res in resSummary){
    ii <- ii + 1
    temp[[ii]] <- data.frame(resSummary[ii])
  }
  openxlsx::write.xlsx(temp, paste(outDataPath, "\\PtsABCRFSummary.xlsx",
                                   sep=""),
                       rowNames=TRUE)
}